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1  Introduction 


This  document  is  the  hnal  report  for  project  FA  9550-05-1-0439,  funded 
by  AFOSR  through  the  DEPSCOR  program.  The  program  manager  was 
Arje  Nachman.  The  research  was  to  investigate  the  dynamics  of  nonlinear 
waves  and  other  flows  in  the  vicinity  of  the  tropopause  using  theoretical  and 
computational  methods.  The  project  was  very  successful,  and  has  resulted 
in  many  new  and  important  results.  The  primary  results  are  documented 
below. 

A  number  of  students  in  various  capacities  are  now  affiliated  with  the 
project.  They  are  listed  below: 


Name 

Degree 

Nick  Jenkins 

MS 

Bob  Arredondo 

PhD 

Surupa  Shaw 

PhD 

Zhexuan  Zhang 

PhD 

Imani  Lugalla 

BS 

Travis  Gline 

BS 

(part-time) 


In  addition.  Dr.  lordanka  Panayotova,  previously  a  postdoc  at  AFRL  at 
Hanscom,  now  on  the  faculty  at  Old  Dominion  University,  is  collaborating 
with  the  PI  on  the  theoretical  work. 

Project  funds  and  matching  funds  were  used  for  two  months  of  salary  for 
the  PI  during  the  summer  for  2005,  2006,  and  2007,  which  expended  all 
funds  originally  allocated  for  PI  salary.  Project  funds  supported  two  months 
of  salary  for  the  consultant,  T.  R.  Akylas,  expending  all  of  this  part  of  the 
original  budget. 

A  one-year  no-cost  extension  was  requested  and  awarded.  This  extension 
of  time  is  necessary  due  mostly  to  timing  difficulties  associated  with  the 
graduate  students. 
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2  Boussinesq  theory  of  waves  interacting  with 
the  tropopause 

Waves  impinging  on  an  idealized  two-layer  Boussinesq  model  of  the  tropopause 
region  are  treated  with  weakly  nonlinear  theory.  The  results  have  appeared 
recently  in  The  Journal  of  Atmospheric  Sciences 
The  waves  are  assumed  to  be  periodic  in  the  horizontal  and  propagate 
with  permanent  form.  Internal  waves  of  permanent  form  have  been  treated 
previously,  most  notably  by  Thorpe  [30]  and  Yih  [36] .  Yih  [36]  showed  that 
the  background  density  prohle  must  be  adjusted  to  account  for  a  nonlinear 
shift  in  the  mean  streamline.  The  correction  to  the  background  density 
results  in  a  second-order  correction  to  the  wavespeed,  analogous  to  Stokesian 
waves  [28].  However,  Yih  [36]  showed  that  the  correction  for  internal  waves 
is  negative^  meaning  that  larger  amplitude  waves  travel  slower,  opposite  to 
that  for  free-surface  waves.  Thorpe  [30]  and  Yih  [36]  both  considered  a 
conhguration  where  the  flow  is  bounded  on  top  by  a  rigid  lid,  resulting  in 
complete  reflection  of  the  internal  waves.  The  value  of  correction  to  the 
wavespeed  depends  strongly  on  the  wave  reflection.  Partial  wave  reflection 
occurs  in  the  lower  layer  of  the  conhguration  considered  below,  also  resulting 
in  an  important  correction  to  the  wavespeed.  The  present  results  show  that 
this  correction  is  also  negative. 

The  waves  are  made  steady  by  choosing  a  coordinate  system  that  moves 
horizontally  with  the  wave  speed  [also  used  by  Stokes  [28],  Long  [17],  and 
Yih  [36]].  The  governing  equations  are  then  reduced  to  Long’s  equation  [17]. 
For  the  case  considered  here  of  Boussinesq  how  with  constant  Brunt- Vaisala 
frequency  (in  each  layer)  and  no  upstream  shear.  Long’s  model  becomes 
linear,  and  is  given  by 

+  =  0,  (1) 

where  5  is  the  vertical  displacement  of  streamlines  from  an  upstream  or 
background  state,  and  k  is  a  constant.  At  hrst  glance,  it  would  seem  that 
the  solution  to  the  two  layer  problem  considered  here  is  merely  the  sinusoidal 
solution 

5  oc  sm.mxsm.nz,  (2) 

where  m  and  n  are  constants  unique  to  each  layer.  This  solution  does  satisfy 
(1)  exactly,  but  only  meets  the  boundary  conditions  at  the  mean  interface. 

^J.  McHugh,  Journal  of  the  Atmospheric  Sciences,  66,  pp.  1845-1855,  2009. 
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An  accurate  nonlinear  solution  must  meet  the  boundary  conditions  at  the 
actual  interfacial  position,  rather  than  the  mean,  which  is  not  known  before¬ 
hand.  As  a  result,  in  a  multilayer  fluid  (2)  is  only  a  linear  solution,  accurate 
for  inhnitesimal  wave  amplitudes  only.  Important  nonlinear  effects  result 
when  the  displacement  of  the  interface  is  included.  Further  nonlinear  effects 
result  when  the  background  density  prohle  is  adjusted  to  match  the  average 
density  prohle  in  the  presence  of  waves. 

It  is  shown  here  that  the  nonlinear  effects  at  the  interface  result  in  higher 
frequency  internal  waves  propagating  throughout  the  huid.  These  waves  are 
harmonic  to  the  incident  wave  only  at  the  interface,  where  they  contribute 
to  steepen  the  wave.  Away  from  the  interface,  the  harmonics  propagate  at 
an  angle  to  the  horizontal  that  is  different  than  the  incident  wave;  the  in¬ 
terface  causes  the  harmonics  to  be  scattered.  For  some  parameter  values, 
the  linear  solution  gives  an  evanescent  wave  in  the  upper  layer,  meaning  the 
vertical  structure  is  not  oscillatory.  For  all  parameter  values,  the  higher  har¬ 
monic  waves  in  both  layers  become  evanescent  when  the  effective  frequency 
is  greater  than  the  Brunt- Vaisala  frequency. 


2.1  Basic  equations 


The  flow  is  assumed  to  be  incompressible,  inviscid,  and  two-dimensional. 
Stratification  is  present  due  to  a  non-diffusing  quantity.  The  flow  is  then 
governed  by  the  Euler  equations,  the  continuity  equation,  and  the  equation 
of  incompressibility.  Long  [17]  reduced  these  equations  to  a  form  now  known 
as  Long’s  model.  Long’s  model  assumes  a  horizontal  reference  flow,  ^0(2:0); 
with  a  density  prohle,  Poi^o),  where  Zq  is  the  vertical  position  in  the  reference 
how.  The  streamlines  may  be  dehected  by  a  disturbance,  often  considered 
to  be  a  barrier  to  the  how,  such  as  a  mountain.  The  derivation  of  Long’s 
equation  is  given  by  Long  [17],  and  will  not  be  repeated  here.  The  resulting 
equation  is 


V2(5  + 


11  dg  ^05 
2q  dz  dz 


(3) 


where 


6  =  z  -  zo, 


(4) 


is  the  streamline  displacement,  and 


9  dpo 

Po  dzo  ’ 


(5) 
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q  =  poul-  (6) 

If  the  Boussinesq  approximation  is  assumed  and  uo  is  taken  to  be  constant, 
then  Long’s  equation  reduces  to 


7^S+^S  =  0, 


Un 


(7) 


where 


1  <9po 
Po  dzo  ■ 


(8) 


The  mean  density  prohle  is  chosen  to  be  continuous,  but  have  a  discon¬ 
tinuous  hrst  derivative,  such  that  the  fluid  exists  in  two  semi-inhnite  layers, 
each  layer  having  a  unique  value  of  the  buoyancy  frequency.  The  kinematic 
interfacial  condition  is  that  the  normal  velocity  must  be  continuous  across 
the  interface.  The  dynamic  interfacial  condition  is  that  the  pressure  must 
be  continuous  across  the  interface.  As  stated  by  Durran  [8],  these  conditions 
can  be  met  by  choosing  5  to  obey 


(5i  —  82  (9) 

and 

^Iz  =  ^2z  (10) 

at  the  interface,  where  (5i  is  the  streamline  displacement  held  in  the  lower 
layer,  and  82  for  the  upper  layer.  A  complete  derivation  of  these  conditions  is 
provided  in  the  appendix.  Note  that  these  are  still  the  fully  nonlinear  inter¬ 
facial  conditions,  and  not  just  the  linear  condition.  Of  course,  the  condition 
must  still  be  met  on  the  actual  interfacial  position,  rather  than  the  mean 
interfacial  position,  a  feature  that  results  in  the  important  nonlinear  effects. 

The  vertical  position  of  the  origin  of  the  coordinate  system  is  chosen  to  be 
at  the  mean  position  of  the  interface,  so  that  the  streamline  that  corresponds 
to  the  interface  is  given  by  ;2o  =  0.  The  displacement  of  the  interface  is  then 
determined  by 

8[x,  z)  =  z.  (11) 

Inverting  this  expression  to  obtain  the  interfacial  shape,  even  for  relatively 
simple  expressions  for  8,  is  nontrivial.  Assume  such  a  solution  to  (11)  is 

2:  =  r][x).  (12) 
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Note  that  (5(a;,  O)  is  not  the  same  as  77(0;),  except  in  the  linear  case.  The 
interfacial  displacement,  77,  is  related  to  the  velocity  held  by  the  familiar 
kinematic  interfacial  condition  for  each  layer: 

r]t  +  -  Sz)r]x  =  So;  (13) 


on  the  interface. 

Expand  (9),  (10),  and  (13)  in  a  Taylor  series  about  the  mean  position  of 
the  interface,  and  insert  into  the  interfacial  conditions  to  obtain 


z=0 


z=0 


V+2^izz 


rf+.  ■■  =  82 


z=0 


+5: 


2z 


z=0 


2zz 


z=0 


77^  H -  (14) 


z=0 


+  5l2 


z=0 


z=0 


V+2^1zzz 


772  +  .  ..  =  §2 


z=0 


+  (^2 


z=0 


z=0 


V+2^2zzz 


77  H - 


z=0 


ht  +  (^1  - 

h^  +  f  1  - 


+  ^Iz 


z=0 


77  H - 


2  =  0 


+  ^2 


2=0 


77  H - 


2=0 


Vx 

Vx  f^2x 


+ 


2=0 


(15) 
h  +  ---  )  (16) 


2  =  0 


+  5; 


2xz 


z=0 


h  +  ---  (17) 


2=0 


2.2  Finite  amplitude  waves 

The  governing  equation  in  each  layer  is  the  linear  Long’s  equation, 

+  ^61  =  0,  (18) 

uf 

V%  +  ^62  =  0,  (19) 

ui 

where  ui  and  U2  are  background  velocities  for  each  layer.  Expand  5  in  a 
power  series  in  e,  the  ratio  of  the  wave  amplitude  to  horizontal  wavelength; 

(5i  =  e(5ii  +  €^612  +  c^(5i3  +  •  •  •  ,  (20) 

and  similar  expressions  for  S2  and  77.  Note  the  subscript  convention;  the  hrst 
index  indicates  the  layer,  while  the  second  indicates  the  order. 

A  coordinate  system  is  chosen  to  be  moving  with  the  wavespeed  (as  yet 
undetermined)  to  make  the  how  steady,  and  the  wavespeed  must  also  be 
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expanded  in  the  wave  amplitude  to  suppress  secular  terms.  Both  features 
are  treated  with 


Ui  —  Cio  [l  +  ecli  -|-  e^ci2  +  ■  ■  ■  ]  . 

(21) 

U2  =  C20  [1  +  CC21  -|-  e^C22  +  ■  ■  ■  ]  . 

(22) 

Ultimately,  Ui  and  U2  must  be  equal  to  the  wavespeed  and  to  each  other,  so 
that  the  flow  is  steady.  However  it  is  convenient  to  maintain  them  separately 
for  now,  and  then  ultimately  equate  them  to  determine  the  wavenumber  in 
the  upper  layer. 

To  expedite  the  correction  to  the  background  density  profile,  f3  is  also 
expanded  in  the  wave  amplitude: 

A  =  Ao  [1  +  e/3ii  +  e^/3i2  ■  ■  ■  ]  5  (23) 

and  a  similar  expression  for  (32- 
The  first  order  governing  equations  are 

V2<5n  +  =  0,  (24) 

*^10 

V252i  +  ^^21  =  0,  (25) 

^20 

and  the  interfacial  conditions  are 

<^11  =  <^21,  (26) 

<^11^  =  ^2iz  (27) 

Vlx  *^21x)  (28) 

on  2:  =  0. 

The  solution  to  (24-28)  is  chosen  to  be  an  upwardly  propagating  incident 
wave  in  the  bottom  layer.  The  linear  interfacial  conditions  require  a  reflected 
in  the  lower  layer  and  a  transmitting  wave  in  the  upper  layer: 


5ii 


He* 


—iniz 


+ 


/  ni  -  n2 
\ni  -h  n2 


^iniz 


/  ni  -  n2 
\ni  +  n2 


^-iniz 


(29) 
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(521  — 


2ni 


^  712 ^—i{mx—n2Z^ 


(30) 

^ni+na/'  ''  ^  ’ 

where  ni  and  n2  are  the  vertical  wavenumbers  in  the  lower  and  upper  layers, 
respectively,  and  A*  is  the  complex  conjugate  of  A.  Note  that  (29)  and  (30) 
are  chosen  to  satisfy  a  radiation  condition  in  each  layer. 

The  corresponding  displacement  of  the  interface  is 


Vi  = 


The  dispersion  relations  are 


2ni 

jii  +  n2_ 

^10  ~ 


^20  ~ 


gPio 

m?  +  nf 

gp2o 

rn?  +  n2 


2  ’ 


2  • 


(31) 

(32) 

(33) 


The  dispersion  relation  for  the  upper  layer  must  give  the  same  wave  speed, 
and  share  the  same  horizontal  wavenumber,  m,  as  the  lower  layer,  which  de¬ 
termines  the  vertical  wavenumber  in  the  upper  layer.  The  hrst  approximation 
to  this  vertical  wavenumber  is 


Tir,  =  m 


P20  _  . 
Ao 


,  2/^20 

'Ao' 


(34) 


2.2.1  Second  order 

The  second-order  governing  equations  are 


+  ^(5i2  =  -2cnV25n 


gf^i 


cfo 

V"(522  +  ^(522  =  -2C2iV^5: 


cfo 

g(^2 


-(5ii, 


21 


no 


-20 


1^21, 


(35) 

(36) 


Before  proceeding  further,  the  correction  to  the  upstream  condition  must 
be  considered  to  determine  Pn  and  /52i-  This  is  achieved  as  described  in  the 
next  section.  The  analysis  shows  that  there  is  no  correction  at  this  order, 
matching  the  conclusions  of  Yih  [36],  and  the  details  will  not  be  given: 


/dll  —  f^21  —  0. 


(37) 


The  only  secular  term  left  in  (24)  and  (25)  is  suppressed  by  choosing  cn  = 
C21  =0,  making  the  second  order  equations  homogeneous. 

The  second-order  interfacial  conditions  are 


(^12  ~  (^22  — 


(^21  r  —  S 


lU 


Vu 


(38) 


Sl2z  —  S22Z  — 


^2lzz  —  ^llzz 


d 


Vl: 


V2x  -  ^12x  -  ^  , 

V2x  -  ^22x  =  ^  (s21zVl  )  , 


on  =  0.  Using  (27)  reduces  (38)  to 


^12  —  ^22  —  0, 


(39) 

(40) 

(41) 


(42) 


on  =  0.  Inserting  (29),  (30),  and  (31)  into  (39)  then  simplifying  gives 


Sl2z  —  ^22z  — 


2ni 

ni  +  n2 


J^2^2mx  ^  j^*2^-i2mx  ^  2AA* 


(43) 


on  =  0.  The  forcing  terms  in  (43)  are  not  resonant,  and  require  a  solution 
of  the  form 

5i2  =  B  sin  {2'mx  +  n^z)  C  sin  (712;^) ,  (44) 

where  B  and  C  still  need  to  be  determined.  The  vertical  wavenumber  in 
(44),  ni2,  is  not  determined  by  the  interfacial  conditions,  and  is  chosen  to 
satisfy  the  governing  equations  in  the  bottom  layer,  reflected  in  the  dispersion 
relation  for  the  bottom  layer: 


^12  =  ■'^1  ~  (45) 

The  term  C*  sin  (712;^)  balances  the  constant  in  (43),  and  represents  a  wave- 
induced  mean  flow,  but  driven  by  the  interfacial  conditions.  The  governing 
equation  determines  the  exponent: 

2  9P10  2,2  f  Aa\ 

7i2  =  =m  +n^.  (46) 

^10 
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Similarly  in  the  upper  layer, 


622  =  -Bsin  {2mx  -  n22z)  +  Csin  (722^), 


^22 


=  Tin  —  3m^, 


722  = 


gp2o 


(47) 

(48) 

(49) 


-20 


where  the  coefficients  have  already  been  chosen  to  satisfy  (42).  Note  that 
the  negative  sign  in  front  of  n22  is  chosen  in  the  upper  layer  to  obtain  an 
upwardly  propagating  harmonic,  meeting  the  radiation  condition  at  the  top. 
Equation  (39)  determines  B  and  C\ 


B 


—Anl 


ni  -  n2 
rti  +  n2 


A' 


C 


2ni 

ni  +  n2 


n\  —  nl 
712  +  722 


2AA*. 


(50) 

(51) 


2.2.2  Correction  to  f3 


The  second-order  solution  must  match  the  background  conditions.  This  is 
achieved  using  the  method  of  Yih  [36] .  The  dehnition  of  6  is  rearranged  to 
obtain 

z  =  zo  +  6(^x,  z^.  (52) 

The  inversion  of  this  equation  is 

(53) 


2:  =  rjyx,p 

An  average  over  one  wave  period  of  (53)  is  the  average  displacement  of  a  line 
of  constant  density  in  the  presence  of  waves.  A  wave  of  permanent  form  must 
have  this  vertical  displacement  equal  to  zero.  If  it  is  not,  then  the  upstream 
density  prohle  must  be  adjusted.  This  adjustment  is  determined  using 


1  dp  1  dp  dp  ^  dp 

Po  dzQ  po  dp  dzo  ^  dzo 


(54) 


In  practice,  (52)  is  inverted  using  the  method  of  successive  approximations 
for  each  layer,  as  in  Yih  [36]  and  Stokes  [28].  The  hnal  result  is 

1321  =  4:AA*nl  ( —  cos  7122:0  -  2  cos  2ni2:o)  ,  (55) 

Vni -hn2/ V7i2  +  722  / 
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P22  =  0.  (56) 

Note  that  the  correction  for  the  upper  layer  is  zero  because  there  is  no  reflec¬ 
tion,  only  the  incident  wave.  The  correction  for  the  lower  layer  is  non-zero 
because  of  the  interaction  between  the  incident  wave  and  its  reflection  from 
the  interface. 


2.2.3  Third  order 


The  third-order  solution  is  pursued  just  far  enough  to  demonstrate  uniform 
validity,  and  thereby  determine  the  second-order  correction  to  the  wavespeed. 
The  third-order  governing  equations  are 


V2<5i3  +  =  -2ci2V2<5n 


c 


ni-, 


10 


91320 


9p22 

—^021- 


c; 
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The  third-order  interfacial  conditions  are 


1 

<^13  ~  <^23  — 

<^212  ~  <^112 

V2  + 

322z  ~  ^122 

’h  +  2 

<^2122  ~  ^1122 

(57) 

(58) 


(59) 


<^13^  —  ^23^  — 


321zz  —  3i1zz 


V2  + 


322zz  —  3i2zz 


1 

'(1  +  2 


^21^^2  ~  <^11222 


vl, 


d 


VSt  +  VSx  ~  '^13x  —  ^  (  <^ll2h2  +  <^122hl  +  )  ’ 


d 

VSt  +  VOx  ~  323x  =  ^  (  <^2l2h2  +  322zVI  +  o<^2l22hl 


(60) 

(61) 

(62) 


on  =  0. 

The  forcing  terms  in  the  interfacial  conditions  are  not  resonant,  nor  will 
they  be  at  any  order,  and  hence  they  do  not  contribute  to  the  second-order 
correction  to  the  wave-speed.  This  is  because  the  wave  speed  is  not  deter¬ 
mined  by  the  interfacial  conditions,  as  interfacial  waves  do  not  exist  when  the 
density  is  continuous.  The  interfacial  forcing  terms  are  still  important  how¬ 
ever,  as  they  result  in  scattered  harmonics  of  the  incident  waves,  discussed 
in  the  next  section. 

The  forcing  terms  in  the  governing  equation  for  the  lower  layer  are  secular, 
as  a  result  of  the  nonzero  value  of  f3i2.  These  terms  are  suppressed,  again 
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following  Yih  [36],  by  multiplication  with  the  expression  for  5i  in  the  lower 
layer,  and  integration  over  a  single  wave  period.  The  resulting  second-order 
correction  to  the  wavespeed  in  the  lower  layer  is 


Ci2  =  2n\AA* 


ni-n2\  f  1 
ni  +n2/  V712  +  722/ 


1  Ani-n2)^\ 

2V  n\  +  nl  ) 


7l2 


712  +  722 


(  1  ^  sin(2ni-7i2)^ 

fan.- 712  J  ir 


/  1  \  sin(2ni  +  712)^' 

V2ni+7i2y  ^ 


(63) 


while  for  the  upper  layer, 

C22  =  0.  (64) 

The  hnal  stage  is  to  determine  n2,  governed  by  the  equality  of  wave-speeds 
in  the  two  layers.  Setting  ui  =  U2  in  (21)  and  (22)  and  keeping  terms  to 
second  order  results  in 


9(^10 
w?  +  n\ 


1  +  e^Ci2 


91^20 
m?  +  n\ 


(65) 


This  expression  along  with  (63)  represents  two  coupled  algebraic  expressions 
for  n2  and  C12,  which  are  determined  numerically. 


2.3  Discussion 

Internal  waves  only  exist  for  frequencies  that  are  less  than  the  buoyancy 
frequency,  N  =  gj3,  as  can  be  shown  from  the  dispersion  relation.  Hence  the 
incident  waves  in  the  lower  layer  must  have  a  frequency  less  than  gj3iQ.  The 
buoyancy  frequency  above  the  interface  may  be  larger  or  smaller  than  below: 
the  buoyancy  frequency  at  the  tropopause  and  the  mesopause  increases  with 
the  vertical  coordinate,  while  the  buoyancy  frequency  at  the  stratopause 
decreases. 

For  the  case  where  N2  <  Ni  (the  stratopause),  it  is  possible  for  the  fre¬ 
quency  of  the  linear  incident  wave  to  be  greater  than  1+2.  The  linear  trans¬ 
mitted  wave  will  be  evanescent  in  the  upper  layer,  with  no  vertical  oscillations 
in  the  upper  layer. 
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For  the  case  where  N2  >  Ni  (the  tropopause  and  mesopause),  the  linear  in¬ 
cident  wave  creates  a  wave  in  the  upper  layer  that  oscillates  with  the  vertical. 
However,  the  harmonics  that  are  created  at  the  interface  may  be  evanescent 
in  either  layer.  The  harmonic  may  be  considered  to  have  an  effective  fre¬ 
quency;  for  example,  the  second  harmonic  in  the  lower  layer  would  have  the 
effective  frequency. 


2m 


a/ (2m)2  -1-  n\ 


iVi. 


(66) 


If  the  effective  frequency  of  a  harmonic  is  greater  than  the  Brunt- Vaisala 
frequency  in  either  layer,  then  the  harmonic  will  be  evanescent  in  that  layer. 
This  transition  occurs  when  the  vertical  wavenumber  becomes  purely  imag¬ 
inary,  for  example,  if  n\2  >  0  in  (45),  then  ni2  is  real  and  the  harmonic  is 
oscillatory.  If  n^2  <  0)  then  the  behavior  is  evanescent. 

Away  from  the  interface,  the  higher  harmonics  with  vertical  oscillation  will 
not  coincide  with  the  primary  mode.  For  example,  the  incident  waves  in  the 
lower  layer  are  traveling  waves  with  wavenumbers  m  and  ni,  as  discussed 
previously.  The  lines  of  constant  phase  for  this  wave  make  an  angle,  ^i,  to 
the  horizontal: 

9i  =  arccos  — ,  (67) 

+  nl 

The  second-order  solution  requires  a  harmonic  with  wavenumbers  2m  and 
ni2,  where  ni2  is  given  by  (45).  This  harmonic  is  a  wave  that  makes  an 
angle,  612-,  with  the  horizontal: 


6*12  =  arccos 


2m 


2m 


a/  {2mY  +  n 


2 

12 


=  arccos 


a/ m^  +  n\  ’ 


(68) 


where  (45)  has  been  used  to  simplify  (68).  Clearly,  the  angle  of  the  inci¬ 
dent  wave  is  different  than  the  angle  of  the  second-order  harmonic  wave;  the 
second-order  harmonic  has  been  scattered  by  the  interface. 

The  third-order  interfacial  conditions  will  result  in  third-order  harmonics 
with  a  higher  effective  frequency  than  the  second-order  harmonic.  Depending 
on  the  parameter  values  of  the  incident  wave,  this  frequency  could  be  greater 
than  the  buoyancy  frequency  in  either  layer,  resulting  in  non-oscillatory  be¬ 
havior,  and  an  evanescent  mode.  Furthermore,  eventually  there  is  a  har¬ 
monic  in  the  Stokes  expansion  that  will  result  in  a  harmonic  frequency  that 
is  greater  than  both  Ni  and  iV2,  and  this  harmonic  will  be  evanescent  in  both 
layers,  as  will  all  higher  harmonics. 
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The  deflection  of  the  interface  due  to  the  second-order  harmonics  has  the 
same  phase  as  the  deflection  of  the  interface  due  to  the  linear  waves  given  by 
(31),  whether  the  harmonic  is  oscillatory  or  not.  The  third-order  harmonic 
will  also  have  this  phase.  As  a  result,  the  nonlinear  wave  at  the  tropopause 
behaves  as  a  Stokes  wave,  where  the  crest  is  sharpened  and  the  trough  broad¬ 
ened. 

However,  the  behavior  of  the  waves  away  from  the  interface  is  significantly 
different.  Only  the  oscillatory  harmonics  extend  away  from  the  interface, 
and  depending  on  the  wavenumbers  for  the  incident  wave,  only  the  first  few 
harmonics  are  likely  to  be  oscillatory.  The  remainder  of  the  harmonics  will 
decay  exponentially  in  each  layer.  Hence  a  horizontal  profile  would  consist 
of  the  sum  of  only  a  few  sinusoidal  components,  quite  different  than  the 
behavior  at  the  interface. 

The  second-order  correction  to  the  wave-speed  and  the  wavenumber  in  the 
upper  layer  are  determined  by  the  two  coupled  nonlinear  algebraic  equations 
given  in  (63)  and  (65).  Values  of  Ci2  and  n2  have  been  determined  numerically 
using  the  bisection  method,  and  are  shown  in  figures  2.3  and  2.3.  Figure 
2.3  shows  that  Ci2  is  negative,  meaning  that  larger  amplitude  waves  travel 
slower  than  infinitesimal  waves.  This  result  was  demonstrated  by  Yih  [36] 
for  internal  waves  bounded  with  rigid  horizontal  surfaces.  The  rigid  surfaces 
in  Yih  [36]  cause  complete  wave  reflection,  which  in  turn  leads  to  a  non¬ 
zero  wave-induced  mean  flow  and  a  displacement  of  isopycnic  lines.  The 
present  results  have  partial  wave  reflection  from  the  interface,  rather  than 
the  complete  reflection  as  in  Yih  [36].  Furthermore,  there  is  an  additional 
contribution  to  the  wave-induced  mean  flow  as  a  result  of  the  interfacial 
conditions,  given  by  the  second  term  in  (44).  This  interfacial  mean  flow 
weakens  the  effect,  but  is  not  significant  enough  to  overcome  the  effect  found 
by  Yih  [36],  and  hence  the  same  mechanism  that  results  in  the  negative 
wavespeed  correction  is  also  at  work  here. 
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3  DNS  of  waves  at  the  tropopause 

The  direct  numerical  simulation  (DNS)  of  waves  interacting  with  an  idealized 
tropopause  was  treated  for  a  variety  of  conditions.  Waves  are  created  beneath 
the  interface  and  allowed  to  impinge  on  the  interface.  The  results  show  that 
a  horizontal  mean  flow  is  created  at  the  interface  that  grows  in  strength 
and  hnally  creates  a  critical  layer.  These  results  have  been  published  in 
Theoretical  and  Computational  Fluid  Dynamics 


3.1  Governing  equations 


The  governing  equations  are  the  anelastic  equations  [16,  19,  34].  The  anelas- 
tic  equations  are  an  approximate  form  of  the  Navier-Stokes  equations  where 
the  effects  of  sound  wave  propagation  have  been  removed.  The  original 
anelastic  equations  are  flawed,  and  the  linear  solution  does  not  match  the 
linear  solution  of  the  fully  compressible  equations.  It  was  later  shown  that 
the  pressure  term  was  incorrect.  Bacmeister  and  Schoeberl  [4]  provide  a  con¬ 
cise  derivation  of  the  corrected  anelastic  equations,  which  have  been  found 
to  accurately  model  atmospheric  dynamics. 

For  large-scale  waves  typically  found  in  the  atmosphere,  viscosity  is  not 
important;  inviscid  flow  is  asssumed  here.  The  inviscid  anelastic  equations 
in  two  dimensions  for  a  perfect  gas  atmosphere  are 


where 


Du  dp* 

(69) 

Dt  dx  ’ 

Dw  dp*  9 

Dt  ~  dz 

(70) 

DO  dO  ^ 

(71) 

~  ’ 

djiu  dpv 
dx  ^  dy 

(72) 

R 

*  a  i  ^  V'’ 

p  =Cp9  {  —  ]  , 

\PoJ 

(73) 

d  d  d 

(74) 

^McHugh,  J.  P.,  Theoretical  and  Computational  Fluid  Dynamics,  22,  ppl07-123,  2008. 
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M,  w  are  the  velocity  components,  x,  2;  are  the  components  of  position,  9  is 
the  potential  temperature,  p  and  9  are  the  base  state  density  and  potential 
temperature,  respectively,  Cp  is  the  specihc  heat  at  constant  pressure,  R  is 
the  gas  constant,  g  is  the  gravitational  constant,  p  is  the  pressure,  and  po 
is  a  constant.  Equations  (69-70)  are  the  momentum  equations,  (71)  is  the 
energy  equation,  and  (72)  is  the  continuity  equation. 


3.2  The  base  state 

The  base  state  is  governed  by  the  perfect  gas  law, 

p  =  pRT, 


(75) 


and  the  equation  of  static  equilibrium, 

dp 
dz 

Note  that  all  base  state  variables  are  denoted  by  an  overline.  An  important 
base-state  parameter  is  the  Brunt- Vaisala  frequency,  dehned  as 


=  -P9- 


(76) 


Ar2  =  1‘P 

9  dz 


(77) 


A  variety  of  base  states  may  be  considered,  including  cases  drawn  from 
observations.  However,  the  base  state  must  satisfy  the  above  base  state 
equations  to  be  valid.  Often  the  base  state  temperature  is  chosen,  and  then 
the  density  and  potential  temperature  are  determined.  The  following  two 
relations,  which  are  derived  directly  from  the  first  law  of  thermodynamics 
for  a  dry  atmosphere  and  the  dehnition  of  potential  temperature  [13],  are 
used  to  determine  p  and  9  once  T  is  chosen: 


^  +  ^  +  4=  =0, 

pTRT 

Oz  T,  9  _Q 

9  T  CpT 


(78) 


(79) 


The  base-state  temperature  for  the  present  results  is  chosen  to  be  constant 

2 

in  each  layer,  making  the  Brunt-Vaisala  frequency  constant  and  equal  to 

Cpl 

with  a  unique  value  in  each  layer.  The  potential  temperature  is  given  by 


9  = 


(80) 
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where  the  exponent  has  a  unique  value  in  each  layer.  The  pressure  must 
be  continuous  at  the  interface  for  equilibrium,  and  density  is  chosen  to  be 
continuous.  The  potential  temperature  is  also  chosen  to  be  continuous  at  the 
interface. 

Another  important  base  state  parameter  is  the  scale  height,  H,  which  is 
the  thickness  of  a  layer  of  gas  of  constant  temperature  over  which  the  density 
changes  by  the  factor,  e.  The  local  value  of  the  scale  height  is  dehned  by 

H  p- 

For  a  constant  temperature  layer,  the  scale  height  is  given  hy  H  =  ^.  The 
above  choices  result  in  a  constant  scale  height  throughout  the  domain. 

3.3  The  computational  approach 

The  velocity,  position,  and  pressure  are  rescaled  using  a  length  scale,  L, 
and  a  velocity  scale,  U.  The  wave  forcing  is  chosen  to  have  a  horizontal 
wavelength.  A,  that  is  equal  to  the  horizontal  length  of  the  domain,  and  this 
lengthscale  is  chosen  as  the  lengthscale  for  rescaling;  L  =  X.  The  flow  is 
initiated  from  rest;  the  initial  mean  flow  is  also  zero.  This  means  that  the 
basic  state  does  not  provide  a  natural  velocity  scale.  Instead,  the  velocity 
scale  is  chosen  to  be  V^.  The  density  and  potential  temperature  are  rescaled 
using  their  respective  base-state  values  at  the  bottom  of  the  domain  (^,  6*o). 
The  dimensionless  equations  are 


Du  dp 

Dt  dx' 

(82) 

Dw  dp  9 

Dt  dz 

(83) 

DO  ^dd  ^ 

+ 

Dt  dz 

(84) 

dpu  dpv 
dx  ^  dy 

(85) 

where  the  circumflex  denotes  a  dimensionless  quantity. 
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The  governing  equations  are  reduced  such  that  the  pressure,  horizontal 
velocity,  and  potential  temperature  are  eliminated  from  the  linear  terms, 
resulting  in 


d  (\ 


V  w;  +  —  -^w 
dz  \h 


+  N^V\w  = 

dt 


d  OAj 
dz  dxi 


B 

-V? 

86) 

where  the  Ai  is  the  sum  of  the  nonlinear  terms  for  the  momentum  equa¬ 
tion, 

Ai  =  Uj^,  (87) 


B  is  the  sum  of  the  nonlinear  terms  for  the  energy  equation. 


and 


R  - 

B  =  Uj  — , 

Xj 

(88) 

(92 

v?=  . 

(95;2 

(89) 

The  boundary  condition  on  the  top  boundary,  z  =  ^,  where  D  is  the  height 
of  the  domain,  is  tc  =  0.  Note  that  a  damping  layer  is  included  near  the  top 
boundary  to  absorb  upward  propagating  waves,  as  will  be  discussed  later. 
On  the  bottom  boundary,  w  is  the  imposed  forcing  velocity,  also  discussed 
later.  The  side  boundaries  are  treated  as  periodic. 

The  order  of  (86)  is  reduced  by  introducing  0: 


t-72  ~  ^ 

oz 


d  dAj 
dz  dxi 


(90) 


Equation  (86)  becomes 


^  +  N^Vlw 

dt  ^ 


(91) 


The  variable  0  is  retained  in  the  calculations,  and  w  and  0  are  determined, 
using  (90)  and  (91),  for  each  time  step.  Note  that  the  nonlinear  terms  that 
contain  a  temporal  derivative  have  been  carefully  chosen  as  part  of  the  def¬ 
inition  of  0.  This  subtle  but  important  step  allows  the  nonlinear  terms  to 
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be  evaluated  by  accurate  well-known  methods,  such  as  Adams-Bashforth,  to 
any  level  of  accuracy. 

Equation  (90)  and  (91),  being  first-order  in  time,  may  be  directly  integrated 
over  the  time  step.  At,  to  obtain 


^  I  -  ~ 

yw  +  —  -^w 


dz  \H 


n+1 


T-r2~  ^  I  -  ~ 

yw  +  —  -^w 


dz  \H 
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fn  +  1 


d  dAi 


Zn  +  l 
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dz  dx 

i 

n  +  l 

B 

dt  =  — 

V? 

Ji^ 

% 

(j)  dt  = 


- 


dt. 


di,  (92) 


(93) 


The  integrals  in  the  linear  terms  in  (92)  and  (93)  are  treated  implicitly  with 
the  Crank-Nicolson  method,  while  the  nonlinear  terms  are  treated  explicitly 
with  the  second-order  Adams-Bashforth  method.  The  result  is 


A 

dz 


n+1 

d  dAj 
dz  dxi 


_2  _  d 
V  w  — 
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0n+l  ^ 
2 


At  = 


f^dA 
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At,  (94) 


In+l 


At  =  -^V? 
9 


3B^  -  B 


n-n 


At.  (95) 


There  is  one  remaining  difficulty;  both  equations  contain  w  and  0  at  the 
leading  time  step,  making  them  coupled.  The  two  equations  may  be  decou¬ 
pled  merely  by  eliminating  algebraically,  resulting  in 
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Equation  (96)  may  be  solved  directly  for  Values  of  may  then  be 

determined  using  either  (94)  or  (95);  (95)  is  used  in  practice  for  efficiency. 
The  continuity  equation  is  used  to  determine  the  horizontal  velocity: 


Ux  = 


Pz  ~ 
Wz  +  —w 

p 


(97) 


This  equation  is  easily  solved  for  once  is  available.  That  leaves 

only  9,  which  can  be  found  directly  from  the  energy  equation, 


nn+l  nn  _ 

e  -0  --- 
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n+l 


+  W' 


At 


-  B 


n—1 


At. 


(98) 


The  continuity  equation  cannot  be  used  to  determine  the  horizontal  average 
of  horizontal  velocity,  {u),  where  the  angle  brackets  indicate  the  horizontal 
average,  since  any  prohle  of  (u)  will  satisfy  continuity.  Instead  the  horizontal 
average  of  the  horizontal  momentum  equation  is  used: 


d{u) 

~W 


1  d 

■z^ipnw). 


(99) 


Note  that  the  above  formulation  is  easily  extended  to  three  dimensions. 

The  spatial  discretization  is  a  spectral  method.  The  horizontal  direction 
is  expanded  in  a  Fonrier  series.  The  vertical  direction  uses  an  expansion  in 
Cardinal  fnnctions,  collocated  on  the  Chebyshev-Gauss-Lobatto  grid.  Details 
of  this  spatial  method  can  be  found  in  Solomonoff  and  Turkel  [27]  and  Boyd 
[6].  Matrix  equations  are  treated  with  direct  methods. 

Note  that  N  is  discontinous  at  the  interface,  which  is  sometimes  incom¬ 
patible  with  a  spectal  method,  leading  to  Gibbs  type  oscillations.  For  the 
present  method,  since  the  spatial  treatment  is  collocation,  the  actual  value  of 
N  at  each  gridpoint  is  used  directly  in  the  discrete  equations.  The  resulting 
solution  for  in  (94)  is  continnous,  even  thongh  N  is  discontinuous,  hence 
Gibbs  oscillations  do  not  appear. 

Waves  are  forced  by  imposing  an  artihcial  vertical  velocity  at  the  bottom  of 
the  domain.  The  imposed  freqnency  and  amplitnde  determine  the  character 
of  the  resulting  wavetrain.  The  forcing  creates  a  wave  packet  with  hnite 
vertical  extent.  The  envelope  function  is  the  raised  cosine.  This  bottom 
forcing  is  given  by 


m)(0,  X,  t) 


sin 


if  i  <  jf 
if  t  >  jf, 


(100) 
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where  a  is  the  forcing  frequency,  a  is  the  forcing  amplitude,  T  is  the  funda¬ 
mental  wave  period,  and  j  is  a  chosen  number  of  fundamental  wave  periods. 

A  damping  layer  is  employed  at  the  top  of  the  domain  as  a  non-reflecting 
lid.  The  damping  layer  is  Rayleigh  friction  [23],  where  the  quantity  j3u  is 
included  in  the  momentum  equations.  The  damping  layer  is  similar  as  that 
used  by  Klemp  and  Lilly  [15],  and  discussed  further  by  Durran  [9].  The 
damping  terms  are  treated  explicitly,  and  any  derivatives  of  (3  that  arise  in 
the  formulation  are  neglected. 

Filtering  is  used  to  avoid  the  accumulation  of  energy  at  the  highest  resolved 
frequencies.  Filtering  is  achieved  here  using  the  sequence  of  spectral  Liters 
discussed  by  Vandeven  [32].  The  filter  is  defined  by 

where  uj  is  the  filter  value,  C,  is  the  frequency,  is  a  dummy  variable  for 
C,  and  p  is  the  integer  order  of  the  Liter.  This  integral  may  be  evaluated 
analytically  to  obtain 

(102) 

This  family  of  Liters  becomes  increasingly  sharp  as  p  increases,  and  were 
proven  by  Vandeven  [32]  to  retain  spectral  accuracy.  In  addition,  some  of 
the  more  common  Liters  can  be  shown  to  be  members  of  this  sequence.  The 
velocity  and  potential  temperature  Lelds  were  Lltered  in  all  directions  at  each 
time  step  for  a  Lxed  value  of  p.  The  value  of  p  =  15  was  generally  found 
to  be  sufficient.  Furthermore,  simulations  of  the  linear  equations  with  and 
without  Lltering  with  this  value  of  p  demonstrated  no  signiLcant  eLects  on 
the  basic  wave. 

All  cases  reported  here  use  a  resolution  of  128  x  128.  A  number  of  cases 
with  higher  resolution,  up  to  512  x  512,  have  been  considered,  and  show  that 
the  results  are  not  signiLcantly  changed.  Similarly,  the  time  step  is  set  to 
At  =  0.1  for  all  results,  after  a  substantial  eLort  to  demonstrate  insensitivity 
to  the  time  step. 

3.4  Results 

Results  are  obtained  for  two  cases;  1)  A  of  10  km  and  2)  A  of  100  km  (A  is  the 
horizontal  wavelength  of  the  forcing).  The  parameter  values  for  case  1  and 
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2  are  shown  in  tables  1  and  2,  respectively.  Case  1  was  chosen  as  typical  of 
internal  waves  in  the  atmosphere,  while  case  2  proved  to  be  more  convenient 
for  studying  large  amplitude  forcing  and  wave  breaking. 

Results  for  case  1  are  shown  in  hgures  2,  3,  and  4.  The  prohle  indicated 
by  the  solid  line  in  hgure  2  is  a  vertical  slice  of  vertical  velocity,  w,  at  three 
time  values.  The  prohle  indicated  by  a  dashed  line  in  hgure  2  is  the  same 
case  repeated  without  nonlinear  ehects.  Note  that  the  horizontal  dashed  line 
shows  the  vertical  position  of  the  interface  in  all  hgures. 

The  wave  packet  that  is  created  by  the  bottom  forcing  is  chosen  to  have 
a  vertical  length  equal  to  two  vertical  wavelengths  of  the  fundamental  wave, 
while  the  distance  between  the  bottom  boundary  and  the  interface  is  ap¬ 
proximately  2.7  times  this  same  vertical  wavelength.  Hence,  the  wave  packet 
for  case  1  ’hts’  between  the  bottom  boundary  and  the  interface.  Figure  2a 
shows  the  wave  packet  is  below  the  interface,  but  beginning  to  interact  with 
it.  The  results  in  hgure  2b  are  for  a  later  time  where  the  wave  packet  has  al¬ 
ready  interacted  with  the  interface,  creating  transmitted  and  rehected  waves, 
both  substantially  weaker  than  the  incident  waves.  In  hgure  2c,  the  rehected 
waves  have  impacted  the  bottom  of  the  domain  and  rehected  again,  however 
the  packet  has  now  lost  much  of  its  coherence. 

Figure  3  shows  the  horizontal  average  of  horizontal  velocity,  (h),  for  the 
same  times  as  hgure  2.  Figure  3a  shows  that  there  is  a  horizontal  mean 
how  aligned  with  the  wave  packet.  The  simulations  show  that  this  mean 
how  moves  with  the  wave  packet,  as  discussed  previously  by  Sutherland  [29] 
using  constant  N  throughout. 

A  diherent  result  is  evident  in  the  average  horizontal  velocity  shown  in  hg¬ 
ure  3b  and  3c;  the  wave  packet  develops  a  horizontal  mean  how,  concentrated 
near  the  interface.  Note  that  the  wave  packet  in  hgure  2b  at  f  =  1000  has 
moved  beyond  the  interface,  yet  the  mean  how  shown  in  hgure  3b  remains. 
Figure  3c  shows  that  this  mean  how  at  f  =  1500  is  largely  unchanged.  Much 
later  times,  when  the  wave  motion  has  completely  dissipated,  show  that  this 
mean  how  remains  permanently. 

These  results  indicate  that  there  are  two  components  to  the  mean  how,  the 
mean  how  that  is  associated  with  the  wave  packet,  moving  with  the  packet, 
and  the  mean  how  at  the  interface  that  becomes  a  permanent  feature  of  the 
how,  even  after  the  waves  are  gone.  The  strength  of  both  components  of 
the  mean  how  depends  strongly  on  the  amplitude  of  the  wave.  However,  the 
character  of  the  mean  how  does  not  change  dramatically  with  wave  ampli¬ 
tude.  For  example,  hgure  5  shows  the  normalized  mean  how  for  case  1  with 
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three  values  of  forcing  amplitude,  a  =  0.01,  0.05,  and  0.1.  The  normalized 
mean  flow  is  dehned  as 

52 

Figure  5  shows  that  this  normalized  mean  flow  is  nearly  identical  for  the 
three  different  amplitudes.  This  is  true  for  both  the  packet  mean  flow,  hgure 
5a,  and  the  interfacial  mean  flow,  hgure  5b. 

The  mean  how  is  obtained  by  solving  (99),  which  is  the  horizontal  average 
of  (69).  The  presence  of  the  interfacial  mean  how  can  be  explained  as  a 
consequence  of  the  change  in  momentum  hux,  (puw).  Integrate  (99)  across 
the  interface  to  obtain 

d  f 

-- Jj>{u)  dz  =  -  (puw)  (104) 

The  momentum  hux,  (puw),  above  the  interface  is  much  small  than  below 
the  interface,  due  to  the  wave  rehection.  The  wave  rehection  is  necessary 
to  meet  the  interfacial  conditions,  even  for  linear  waves,  as  discussed  in  the 
previous  section.  The  decrease  in  (puw)  across  the  interface  makes  the  right- 
hand-side  of  (104)  positive,  which  can  only  be  balanced  by  an  increase  in  (h); 
a  horizontal  acceleration.  More  generally,  since  the  side  boundary  conditions 
are  periodic,  and  the  how  inviscid,  then  there  is  no  net  force  to  counter  the 
change  in  momentum  caused  by  the  wave  rehection.  The  horizontal  acceler¬ 
ation  appears  instead.  Further  evidence  is  provided  by  the  normalized  mean 
velocity,  hgure  5,  which  shows  that  the  strength  of  the  mean  how  depends 
approximately  on  the  square  of  the  amplitude,  as  does  the  momentum  hux. 

The  behavior  in  the  region  beneath  the  interface  changes  dramatically  dur¬ 
ing  the  evolution  of  the  wave  packet,  as  can  be  seen  in  hgure  2.  In  hgure 
2a,  the  wave  packet  is  well-dehned,  while  in  hgure  2c,  the  motion  is  reduced 
to  a  standing  wave  oscillation.  This  feature  is  more  clearly  shown  in  hgure 
4  with  contours  of  the  vertical  velocity  at  times  that  correspond  to  hgures  2 
and  3.  Figure  4a  shows  the  coherent  wave  packet  in  the  lower  layer,  while 
hgure  4c  shows  that  the  wave  packet  is  now  a  standing  wave,  with  the  vertical 
wavelength  equal  to  the  vertical  distance  from  the  bottom  to  the  interface. 

Once  the  standing  wave  oscillation  begins  in  earnest,  the  strength  decreases 
as  wave  energy  is  transmitted  through  the  interface.  The  resnlting  behavior 
in  the  npper  layer  dnring  this  period  is  a  temporally  decaying  mode  that  is 
evanescent  in  the  npper  layer,  meaning  that  the  wave  amplitude  decreases 
with  increasing  altitnde.  Note  that  an  evanescent  mode  does  not  exist  if 
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the  waves  are  allowed  to  fill  the  entire  domain  and  become  fully-developed, 
as  indicated  in  the  previous  section.  Evanescent  behavior  can  be  seen  at 
earlier  times,  before  the  standing  oscillation  is  predominant.  For  example, 
the  behavior  in  the  upper  layer  in  Figure  2b  shows  that  the  transmitted  wave 
has  an  evanescent  profile. 

Contours  of  total  potential  temperature  (not  shown),  which  are  often  used 
to  indicate  wave  overturning  [33],  clearly  show  that  the  waves  are  not  over¬ 
turning  for  case  1  for  the  amplitudes  mentioned  above.  Higher  amplitude 
forcing  has  been  found  to  result  in  a  shear  flow  instability  at  the  bottom  of 
the  domain,  as  a  result  of  the  mean  flow  created  near  the  wavemaker.  This 
phenomenon  makes  it  difficult  to  create  large  amplitude  waves  for  the  param¬ 
eter  values  of  case  1.  Hence  the  results  of  case  1  correspond  to  a  relatively 
small  amplitude.  This  can  be  seen  in  figure  2,  which  shows  profiles  of  w  for 
linear  (dashed  line)  as  well  as  nonlinear  cases.  The  difference  between  linear 
and  nonlinear  is  largest  in  figure  2b  below  the  interface,  where  incident  waves 
and  reflected  waves  are  interacting,  but  overall  the  differences  are  not  large. 

The  results  for  case  2  are  shown  in  figures  6-17  for  three  forcing  amplitudes. 
Figure  6,  7,  and  8  show  profiles  of  vertical  velocity,  mean  horizontal  velocity, 
and  contours  of  vertical  velocity,  respectively,  for  three  time  values  and  with 
a  forcing  amplitude  of  0.1  (as  for  case  1).  Figure  9  shows  contours  of  total 
potential  temperature  for  the  same  times.  The  interface  is  again  indicated 
with  the  dashed  line. 

The  times  for  figures  6a  and  6b  correspond  to  when  the  wave  packet  is 
interacting  with  the  interface.  The  vertical  position  of  the  wave  packet  can 
be  seen  clearly  in  the  contours  of  vertical  velocity,  shown  in  figure  8.  This 
case  (case  2)  does  not  leave  behind  a  strong  standing  wave,  as  did  case  1.  The 
transmitted  wave  packet  remains  coherent  and  continues  to  ascend,  growing 
in  amplitude.  Finally,  convective  overturning  occurs,  as  can  be  seen  in  the 
contours  of  total  potential  temperature  shown  in  figure  9c.  The  waves  ’break’ 
at  this  point  in  the  simulation.  With  the  chosen  filter  (p  =  15),  the  simulation 
cannot  continue.  Stronger  filtering  would  be  required. 

Note  in  figure  7c  that  the  mean  flow  at  the  interfacial  altitude  has  not 
changed  from  the  mean  flow  in  figure  7b  at  an  earlier  time.  The  wave  packet 
has  moved  beyond  the  interface  and  left  this  interfacial  mean  flow  behind,  as 
with  the  previous  case.  The  strong  mean  flow  above  the  interface  in  figure 
7c  is  moving  with  the  wave  packet. 

The  contours  of  total  potential  temperature  for  the  earlier  times,  shown 
in  figure  9a  and  9b,  show  that  the  wave  steepens  considerably  during  the 
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interaction  with  the  interface.  Figure  9b  shows  in  fact  that  the  wave  nearly 
overturns;  closer  examination  shows  that  it  does  not.  This  wave  steepening 
appears  to  be  directly  related  to  the  presence  of  the  mean  flow,  as  will  be 
seen. 

Results  for  a  higher  amplitude,  a  =  0.15,  are  shown  in  hgures  10-13.  The 
evolution  of  the  wave  packet  is  very  similar  with  this  larger  amplitude.  In 
particular,  hgure  11  shows  that  the  interfacial  mean  flow  has  the  same  form 
as  in  hgure  7,  except  now  is  stronger.  A  signihcant  difference  is  evident  in 
the  contours  of  total  potential  temperature,  shown  in  hgure  13.  Figure  13a 
shows  that  the  wave  is  overturning  below  the  interface.  Figure  13b  indicates 
that  breaking  has  occurred  below  the  interface,  and  the  hltering  has  managed 
to  suppress  the  small  scale  activity  that  is  created  during  the  breaking  event. 

The  wave  packet  has  moved  past  the  interface  in  hgure  13c,  and  has  reach 
an  amplitude  where  breaking  is  occurring  again.  The  position  of  the  wave 
packet  is  evident  in  the  prohles  of  w  in  hgure  10,  and  contours  of  w  in  hgure 
12.  Note  that  the  altitude  for  breaking  is  lower  than  the  altitude  for  breaking 
in  hgure  9c  for  the  case  with  d  =  0.1. 

Results  for  an  even  higher  amplitude,  a  =  0.2,  are  shown  in  hgures  14- 
17.  Once  again,  the  only  diherence  is  the  value  of  the  forcing  amplitude;  all 
other  parameters  are  set  to  the  values  for  case  2  listed  in  table  2.  The  wave 
dynamics  are  again  similar,  except  now  the  wave  breaks  more  violently  below 
the  interface.  This  can  be  seen  in  the  contours  of  total  potential  temperature, 
hgure  17.  The  hltering  is  not  sufficient  to  control  this  breaking  event,  and 
the  smallest  scales  are  overwhelmed.  The  simulation  cannot  continue  without 
added  hltering. 

It  can  be  seen  in  hgures  14-17  that  the  wave  packet  has  been  transmitted, 
and  is  developing  another  breaking  event.  This  is  most  evident  in  hgure  17b, 
which  shows  the  contours  of  total  potential  temperature  as  they  steepen. 

The  wave  breaking  below  the  interface  can  be  explained  as  a  consequence 
of  the  interfacial  mean  how.  The  horizontal  phase  velocity  of  a  linear  wave 
in  this  rescaled  system,  neglecting  the  inverse  of  the  scale  height,  is 

Cp  =  a,  (105) 

It  is  well-known  that  an  internal  wave  attempting  to  transit  a  mean  shear 
how  will  experience  a  critical  layer  at  the  place  where  the  parallel  component 
of  wave  speed  is  equal  to  the  local  mean  velocity.  The  forcing  frequency,  and 
therefore  the  horizontal  wavespeed,  for  case  2  is  approximately  unity.  The 
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interfacial  mean  flow  with  ci  =  0.1  has  a  maximum  value  of  approximately 
unity,  as  shown  in  hgure  7b.  This  mean  flow  apparently  caused  the  wave 
steepening  shown  in  hgure  9b,  but  was  not  strong  enough  to  cause  overturn¬ 
ing. 

The  maximum  value  of  the  interfacial  mean  how  with  d  =  0.15  was  greater 
than  unity,  and  greater  than  the  horizontal  wave  speed.  This  mean  how  is 
strong  enough  to  overturn  the  waves,  as  is  evident  in  hgure  13.  And  hnally, 
with  a  =  0.2,  the  mean  how  is  much  greater  than  unity,  resulting  in  a  much 
more  violent  breaking  event. 

One  hnal  comment  concerns  the  location  of  breaking.  All  simulations  to 
date  show  that  the  interfacial  mean  how  causes  breaking  below  the  interface, 
never  above.  In  fact,  the  wave  amplitude  immediately  above  the  interface 
is  greatly  reduced  by  interfacial  rehection  or  wave  breaking.  This  indicates 
that  wave-induced  turbulence  in  the  atmosphere  is  expected  in  the  region 
below  the  tropopause,  not  above.  This  is  in  contrast  to  the  conclusions  of 
Van  Zandt  and  Fritts  [31],  who  argued  that  wave  breaking  and  turbulence 
are  expected  above  the  tropopause,  rather  than  below. 


(a)  i  =  500  (b)  i  =  1000  (c)  i  =  1500 

Figure  3:  Prohles  of  vertical  velocity,  w  =  for  case  1  with  d  =  0.1 

(the  dashed  prohle  is  linear) 


(a)  i  =  500  (b)  i  =  1000  (c)  i  =  1500 

Figure  4:  Mean  horizontal  velocity,  {«),  for  case  1  with  ci  =  0.1 


Vertical  velocity 


Vertical  velocity 


Vertlcol  velocity 


(a)  i  =  240  (b)  i  =  270  (c)  i  =  340 

Figure  5:  Prohles  of  vertical  velocity,  w  =  -^=,  for  case  2  with  ci  =  0.1 

(dashed  prohle  is  linear) 
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(a)  i  =  240  (b)  i  =  270  (c)  i  =  340 

Figure  6:  Mean  horizontal  velocity,  (u),  for  case  2  with  a  =  0.1 
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4  Non-Boussinesq  internal  waves 

Internal  waves  in  a  non-Boussinesq  atmosphere  have  an  exponentially  grow¬ 
ing  amplitude,  becoming  unbounded  in  an  unbounded  domain.  Weakly 
nonlinear  non-Boussinesq  theory  must  contend  with  this  unbounded  linear 
behavior,  and  result  in  higher  harmonics  that  grow  exponentially  with  an 
even  faster  rate.  This  problem  is  an  impediment  to  further  theory  of  non- 
Boussinesq  waves.  An  alternative  expansion  has  been  found,  and  has  been 
published  in  SIAM  Journal  of  Applied  Math 

Internal  waves  of  permanent  form  are  horizontally  periodic  waves  where 
the  wave  amplitude  does  not  evolve.  Weakly  nonlinear  theory  of  internal 
waves  of  permanent  form  in  a  continuously  stratihed  fluid  was  previously 
considered  by  Thorpe  [30]  and  Yih  [36],  both  treating  incompressible  inviscid 
flow  between  two  rigid  horizontal  surfaces  assuming  constant  Brunt- Vaisala 
frequency.  Thorpe  [30]  showed  that  the  waves  generate  a  mean  flow,  but 
he  did  not  determine  a  nonlinear  contribution  to  the  wave  speed.  Yih  [36] 
determined  the  wave-induced  mean  flow  by  adjusting  the  background  state, 
and  found  a  second-order  correction  to  the  wavespeed. 

Yih’s  results  are  surprising,  and  show  that  the  second-order  wavespeed 
correction  is  negative.  A  negative  value  means  that  larger  amplitude  waves 
travel  slower  than  small  amplitude  waves,  opposite  the  case  for  most  wave 
systems,  such  as  free  surface  waves.  Yih  only  reported  results  with  a  hori¬ 
zontal  wave-number  of  unity. 

Both  Thorpe  and  Yih  used  a  straight-forward  expansion  in  the  wave  am¬ 
plitude;  for  example,  Yih  [36]  expanded  the  streamfunction,  tjj,  using 

^jJ  =  ^jJo  +  e^jJl  +  e^'ip2  ^ -  (106) 


where  e  <<  1  is  the  wave  amplitude.  The  linear  solution  has  the  form 


ipi  =62^  sin  kx  sin  mz. 


(107) 


where  k  and  m  are  the  horizontal  and  vertical  wavenumbers,  respectively,  and 
f3  is  the  Boussinesq  parameter.  The  product  of  the  constant  wave  amplitude 
and  the  exponential  growth,  eea^,  is  often  thought  of  as  an  exponentially 
growing  wave  amplitude.  This  characteristic  is  well-known  in  atmospheric 
waves  [13],  and  all  non-Boussinesq  internal  waves. 

^J.  McHugh,  SIAM  Journal  of  Applied  Math,  66,  pp.  1845-1855,  2009. 
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The  second-order  solution  must  contend  with  nonlinear  quadratic  terms, 
and  Yih’s  resulting  expression  for  'ip2  contains  a  term  similar  to 

ip2  =  ■  ■  ■  sin  2kx  sin  2mz  •  •  •  ,  (108) 


along  with  many  other  terms.  Note  that  the  exponential  growth  with  altitude 
of  1^2  is  twice  the  exponential  growth  of  tjji.  The  third  order  equation  contains 
cubic  terms,  resulting  in  the  behavior  and  hnally  the  order  solution 

77.  [3 

would  contain  behavior.  Clearly  the  rate  of  growth  is  increasing  with 
each  term.  This  straight-forward  expansion  in  wave  amplitude  does  not 
converge  at  any  altitude,  but  is  asymptotic  if  the  top  boundary  is  a  rigid 
lid.  If  different  top  boundaries  are  treated,  the  results  with  this  expansion 
are  not  valid. 

An  alternative  expansion  that  is  more  generally  valid  is  considered  here. 
The  new  expansion  is  suggested  by  the  following  ’model’  equation,  obtained 
by  arbitrarily  discarding  terms  from  the  governing  equation: 

Uz  —  u  +  =  0,  (109) 


where  u 


u{z).  The  exact  solution  is  readily  obtained: 

ee^ 


u  = 


1  -|-  ee"”  ’ 


(110) 


where  e  is  a  constant. 

Let  (  =  ee^  for  convenience,  and  the  exact  solution  becomes 

“  =  TTc' 

Note  the  singular  point  at  C  =  ~L  A  McLauren  series  expansion  of  (111)  is 

CXD 

u  =  J2{-  =  ee"  -  +  ■  ■  •  .  (112) 

n=l 


The  presence  of  the  singular  point  limits  the  validity  of  (112)  to  C  <  1-  This 
McLauren  expansion  is  analogous  to  the  expansion  used  by  both  Thorpe  [301 
and  Yih  [36]. 

An  expansion  without  this  limitation  is  a  Laurent  expansion  about  C  =  —  1 
of  the  form 


OO 


n=0 


1 


(113) 
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The  coefficients  can  be  found  by  merely  rearranging  the  exact  solution: 


u  = 


1  +  C 


=  1 


1  +  C’ 


(114) 


resulting  in  oq  =  1,  Oi  =  —  1,  and  =  0  for  n  >  1.  Hence  the  exact  solution 
to  the  model  equation  is  equivalent  to  this  Laurent  expansion  about  C  =  —  1- 
The  results  of  this  model  equation  suggest  that  a  Laurent-type  expansion 
may  provide  a  solution  that  is  more  generally  valid,  and  such  a  solution  is 
considered  here.  The  dependent  variable  will  be  expanded  in  a  series  of  the 
form 


OO 


E 


ee 


1  -|- 


p 


(115) 


where  a  is  a  constant  yet  to  be  determined,  and  p  is  an  integer.  For  the 
nonlinear  case,  only  the  first  few  terms  of  this  expansion  can  be  determined, 
due  to  the  complexity  of  the  equations.  For  this  reason,  the  results  are  only 
valid  for  small  amplitude,  as  with  Yih’s  previous  theory. 

The  results  found  with  the  new  expansion  show  that  the  second-order  cor¬ 
rection  is  positive  or  negative,  depending  on  wave  parameters.  Furthermore, 
the  wave  amplitude  with  this  new  expansion  no  longer  shows  an  unbounded 
exponential  growth  with  altitude.  Instead  the  wave  amplitude  experiences 
an  exponential  growth  with  altitude,  but  then  asymptotically  approaches  a 
constant  value.  This  behavior  of  wave  amplitude  is  often  called  saturation 
[11,  25,  26],  and  waves  whose  amplitude  no  longer  grows  with  altitude  are 
saturated  waves.  The  results  with  the  new  expansion  show  the  same  general 
trends  as  Yih’s  original  expansion. 

This  new  result  is  discussed  below  in  section  5.  Section  2  discusses  the 
equations  governing  the  problem,  section  3  gives  the  linear  solution,  showing 
that  it  is  identical  to  the  well-known  linear  solution,  and  section  4  gives  some 
of  the  details  of  the  weakly  nonlinear  analysis. 


4.1  Basic  equations  and  Long’s  equation 

The  flow  is  assumed  to  be  incompressible,  inviscid,  and  two-dimensional, 
governed  by  the  Euler  equations,  the  continuity  equation,  and  the  equation 
of  incompressibility: 

Du  dp 

^  Dt  dx' 
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(117) 


Dw 


dp 


^  Dt  dz 
du  dw 
dx  dz 

^  =  0 
Dt 


(118) 

(119) 


where  {u,  w)  are  the  velocities  in  the  {x,  z)  directions,  respectively,  p  is  the 
pressure,  p  is  the  density,  and  g  is  the  gravitational  constant. 

If  a  coordinate  system  is  chosen  to  move  with  the  waves,  then  the  flow  is 
steady,  and  Longs  equation  (3)  may  be  used.  The  system  is  made  dimension¬ 
less  using  a  velocity  scale,  U,  to  be  dehned  later,  a  constant  value  of  density, 
Poo,  and  the  layer  thickness,  d: 


^  d  ^  X  ^  z  ^  p 

5  =  -,  x  =  z  =  p=  — . 

CL  Cl  d  Pqq 


Dropping  the  circumflex.  Long’s  equation  becomes 


r  1 1  dg 


d5 

2qdzo\^~^z~ 
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F^ul 


(5  =  0, 


where 


f^2L 

VP’ 


is  the  Fronde  number. 
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4.2  Linear  solution 


Equation  (3)  requires  the  flow  to  be  steady.  Steady  flow  with  traveling  waves 
is  achieved  by  choosing  the  coordinate  system  to  be  moving  with  the  hori¬ 
zontal  wavespeed,  c.  This  approach  is  adopted  for  the  linear  case  by  choosing 
Fj-Uq  =  c.  Note  that  the  nonlinear  theory  will  include  a  wave-induced  mean 
flow,  and  Uq  must  account  for  part  of  this,  as  will  be  seen. 

Linearize  Long’s  equation  to  obtain 


Expand  6  with 
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where  are  yet  to  be  determined.  Insert  (124)  into  (123),  and  collect  on 
powers  of  the  quantity  A  successful  choice  for  a  is 


2‘ 


The  resulting  coefficient  of  is 


+  {p  —  l)/50P2  + 


^+p{p-2)^ 


{p  -  ^)P<Pp-iz  +  (P  -  l)  (2p  -  3)  y  0P-1 


(125) 


(126) 


for  all  p  >  1. 

A  solution  for  all  p  is 


<pp  =  sin  kx  sinmz^ 

where  k  and  m  are  constants.  Equation  (126)  at  each  order  results  in 


(127) 
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P 

+  m2  +  f  ’ 


(128) 


the  familiar  dispersion  relation. 
The  hnal  result  is 


5  = 


•  oo 


66  2' 


^  I  1  +  662^ 


sin  kx  sin  m2;, 


(129) 


where  cos  mz  may  be  substituted  for  sin  m2:.  This  solution  may  be  verihed 
by  direct  substitution. 

Each  p  >  1  could  also  include  a  homogeneous  solution,  with  two  arbitrary 
constants.  These  homogeneous  solutions  have  the  form 

(pp  =  sin  kx  sin  mz.  (130) 


Note  that  these  homogeneous  solutions  beyond  p  =  1  have  exponential  decay 
with  altitude,  hence  these  homogeneous  components  do  not  compromise  the 
uniform  validity  of  the  result  of  the  next  section. 
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The  solution  in  (129)  is  identical  to  the  traditional  linear  solution  with  an 
exponential  growing  wave  amplitude,  as  can  be  demonstrating  by  quoting 
the  well  know  expansion, 


e 


a 


(131) 


4.3  Weakly  nonlinear  waves  with  rigid  boundaries 

A  continuously  stratihed  layer  of  fluid  with  rigid  top  and  bottom  will  now  be 
treated  with  the  expansion  of  the  previous  section.  This  is  the  same  problem 
treated  by  Yih  [36]  and  Thorpe  [30] . 

The  governing  equation  is  (3).  Dehne  Q: 
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ee 


1  +  ee° 


(132) 


Expand  5  in  a  power  series: 

S  =  Q(pi  +  Q‘^(p2  +  ■  ■  ■  )  (133) 

where  0^  are  not  the  same  as  the  linear  solution.  Note  that  the  analysis  is 
aided  by  the  relation 
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(134) 


Two  features  are  anticipated;  1)  the  waves  will  create  a  mean  flow  at  second 
order,  and  2)  the  wavespeed  must  be  adjusted  at  second  order  to  obtain  a 
uniformly  valid  solution.  In  the  linear  theory  of  the  previous  section,  the 
upstream  velocity,  Uq,  was  chosen  to  be  the  constant  wavespeed,  c,  making 
the  waves  of  permanent  form  steady.  For  the  present  nonlinear  theory,  a 
coordinate  system  is  again  chosen  to  move  with  the  waves,  however  choosing 
uo  to  be  constant  does  not  provide  a  success  result.  Instead,  the  wavespeed 
and  part  of  the  wave-generated  mean  flow  are  merged  together  as  Uq,  and 
then  this  is  expanded  in  the  same  manner: 


FrtlQ  Cg 


1  -|-  QCi  -\-  Q^C2  +  ■  ■  ■ 


(135) 
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where  the  c^-’s  are  constants.  The  constant,  Cq,  will  be  equal  to  the  linear 
wavespeed,  and  the  remaining  c/s  will  be  chosen  to  obtain  a  uniformly  valid 
solution. 

The  work  of  Yih  shows  that  the  upstream  density  prohle  must  be  adjusted 
so  the  the  density  prohle  in  the  presence  of  waves  results  in  the  desired  prohle. 
This  correction  is  included  at  this  early  stage  by  expanding  (3: 


Po 


1  +  Q(3i  +  Q‘^P2  +  ■  ■  ■ 


(136) 


where  the  /d^-’s  also  constants. 


4.3.1  First  order 

Insert  (133),  (135),  and  (136)  into  (3)  and  collect  the  coefficient  of  Q  to 
obtain 
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If 


(138) 


matching  the  value  of  the  linear  solution,  then  (137)  becomes 


1 

^1^ 

01  =  0. 

(139) 

The  solution  to  (139)  is 

01  = 

sin  kx  smmz, 

(140) 

where  k  is  constant. 

m  =  nn 

(141) 

for  any  integer,  n,  and  (128)  determines  Cq.  Note  that  (139)  along  with  (141) 
satishes  the  boundary  conditions  at  2:  =  0, 1. 


37 


4.3.2  Second  order 

Using  the  choice  for  a  given  by  (138),  the  second-order  equation  is 

V^02  +  f3o4>2z  +  ^02  =  Po<Plz  +  (142) 

Cq  4 

4>il  +  (^4>iz  + 

—Cl  ~  2^01 

^  *^0 

+/5l  /3o0l2  +  ~  ^01  • 

^  *^0 

Note  that  the  hrst  two  terms  on  the  right-hand-side  are  part  of  the  linear 
solution. 

The  value  for  Pi  is  chosen  to  match  the  upstream  conditions,  and  has  been 
found  to  be  zero,  in  agreement  with  the  results  of  Yih  [36] .  Details  of  this 
correction  are  delayed  until  later,  when  the  correction  is  not  zero.  The  value 
of  Cl  is  chosen  to  assure  a  uniformly  valid  solution,  and  has  also  been  found 
to  be  zero,  in  agreement  with  both  Thorpe  [30]  and  Yih  [36] . 

A  homogeneous  solution  to  (142)  must  be  included  in  order  to  meet  the 
boundary  conditions.  The  necessary  solutions  are 


^0  2  r  7  • 

e  ‘2  a  cos  jiz  +  b  sm  jzz  , 

(143) 

and 

a  cos  7;^ -|- 6  sin  72;  sm.2kx. 

(144) 

where 

H  —  Lro - 

(145) 

V  =  Go  -  4V  - 

(146) 

If  7^  <  0  or  <  0,  then  a  hyperbolic  form  is  chosen  for  (143)  or  (144), 
respectively. 
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4.3.3  Third  order 

The  third-order  equation  is 

-|-  2/3o03z  +  ^  ~  ‘^Po4>2z  +  6^02  —  (147) 

—  00  01x02a;  +  ^0U  +  ^^^22  +  ^02  ” 

—  2C2  00012  +  ^01  ~  ^01 

^  *20 

+02  00012  +  ^01  ~  ^01  ; 

^  *20 
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This  equation  will  be  used  to  determine  C2;  no  attempt  will  be  made  to 
determine  03  in  its  entirety. 

Before  finding  C2,  it  is  necessary  to  determine  02  so  that  the  upstream 
density  prohle  matches  the  density  prohle  in  the  presence  of  waves  averaged 
over  one  wavelength.  Following  Yih  [36],  the  dehnition  of  6  in  (4)  and  the 
above  solution  to  second  order  are  used  to  obtain 


Z  —  Zq  Q(j)l  +  (148) 

This  is  a  nonlinear  algebraic  equation  for  the  shape  of  a  streamline,  2:  = 
ri{x,Zo),  for  a  chosen  value  of  Zq.  A  streamline  in  this  scenario  is  also  a 
line  of  constant  density,  and  the  relationship  between  p  and  2  upstream 
could  be  used  to  eliminate  Zq,  resulting  in  2  =  ri{x,p).  The  inversion  of  this 
relationship  is  the  density  in  the  disturbed  held,  p{x,z). 

Equation  (148)  is  inverted  using  the  method  of  successive  approximations, 
as  in  Stokes  [28]  and  Yih  [36] ,  and  then  averaged  over  a  wavelength  to  obtain 


V  =  zo  + 
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(^hi  +  hs  j  e  ^2“  ^0  (.Qg 


+  (^hi  +  F3{p,)e  smp,zo  —  h3cos2mzo  —  h4sm2mzo  —  hi\,  (149) 


where  expressions  for  the  hj  and  F3  are  provided  in  Appendix  II. 
The  correction  is  now  found  using  the  approach  of  Yih; 

1  dpo  1  dp  dfj  ^  drj 

Po  dzQ  p  dfj  dzQ  dzo  ’ 


(150) 


which  may  be  further  evaluated  by  taking  a  derivative  of  (149).  The  effect 
is  only  felt  at  second  order,  and  is  accounted  for  by  choosing  02: 


02  —  00^ 


/iF3(/i)  +  yj  cos/i2  +  (^yF3(/i)  - /ij  sin/iz  (^hi  +  hs'^e  ^2 


'00 


00^3  +  2mh4  ]  cos2mz  +  (  2m/i3  —  00^4  ]  sin2m2  —  0ohi  >.  (151) 


The  last  stage  of  the  second-order  solution  is  to  determine  C2,  found  by 
expanding  the  right-hand-side  of  (147)  and  suppressing  secular  terms.  Equa¬ 
tion  (147)  is  multiplied  by  smmz  and  integrated  over  the  vertical  domain. 
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The  integral  of  all  terms  involving  03  is  exactly  zero,  as  can  be  demonstrated 
using  Green’s  second  identity  and  integration-by-parts,  making  sin  mz 
the  correct  factor.  The  hnal  result  is 


T2  -|-  4T3  —  PqTi 

n 


(152) 


where  Ti,  T2,  T3,  and  T4  are  provided  in  appendices.  Note  that  Ti  results 
from  the  nonlinear  terms  in  (147),  T2  from  the  terms  involving  02,  T3  from 
the  20o02z  +  6^02  — ^01  in  (147),  and  T4  from  the  terms  involving  C2. 


4.4  Results 

The  complicated  expression  for  C2  given  by  (152)  has  been  evaluated  numer¬ 
ically  for  a  wide  variety  of  parameters.  Results  for  two  modes,  n  =  1,2,  are 
shown  in  hgures  4.4  and  4.4  using  a  solid  line  with  0o  =  0.1.  Note  that  the 
dimensionless  value  of  0o  for  atmospheric  flows  is  typically  between  0.01  and 
0.1,  while  flows  in  geophysical  bodies  of  water  may  have  a  wider  variety  of 
values. 

Figure  4.4  shows  that  for  k  ^  1,  C2  is  negative.  This  agrees  with  the  previ¬ 
ous  results  of  Yih  [36],  who  only  reports  results  for  k  near  unity.  For  larger 
k,  hgure  4.4  shows  that  C2  is  positive.  Figure  4.4  shows  that  C2  for  the  second 
mode,  n  =  2,  is  positive,  except  for  very  long  waves.  The  physical  conse¬ 
quence  of  a  positive  C2  is  that  higher  amplitude  internal  waves  travel  faster, 
often  called  ’amplitude  dispersion’.  Amplitude  dispersion  exists  in  many 
other  waves  systems,  in  particular,  free-surface  waves.  Amplitude  dispersion 
is  important  for  the  existence  of  certain  types  of  solitary  waves.  Previous 
studies  of  solitary  waves  in  a  continuously  stratihed  flow  between  horizon¬ 
tal  walls  was  considered  by  Benjamin  [5].  Benjamin  [5]  demonstrated  that 
solitary  waves  exist  for  this  conhguration,  but  he  assumed  the  existence  of 
amplitude  dispersion  (pointed  out  by  Yih  [36]).  It  was  Yih’s  work  showing 
a  negative  C2  that  was  thought  to  eliminate  this  possibility.  However,  the 
present  results  indicate  that  amplitude  dispersion  does  exist,  and  solitary 
waves  may  form  in  this  conhguration.  Solitary  waves  can  form  in  stratihed 
how  under  other  conditions,  as  reviewed  recently  by  Pelinovsky,  et.  ah  [22]. 

A  direct  comparison  with  the  straight-forward  amplitude  expansion  of  Yih 
may  be  obtained  by  expanding  Q,  as  dehned  in  (6),  using 

Q  =  ■  ■  ■  (153) 
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and  deleting  all  but  the  first  term.  The  expansion  in  Q  is  then  identical  to 
that  of  Yih,  and  the  present  results  reduce  to  exactly  the  formula’s  of  Yih. 
The  value  of  C2  for  Yih’s  expansion  are  obtained  from  (152)  by  deleting  the 
quantity,  T3.  The  results  without  T3  are  also  plotted  in  hgures  4.4  and  4.4 
with  a  dashed  line,  and  show  the  same  trends  as  the  present  results.  For  all 
cases  considered,  the  present  results  predict  a  wave-speed  correction  that  is 
greater  than  that  predicted  by  Yih’s  expansion. 

Figures  4.4  and  4.4  show  that  C2  becomes  unbounded  for  isolated  values  of 
k.  These  values  correspond  to  7  =  ±n7r,  where  7  is  dehned  in  (146)  as  part 
of  the  second-order  homogeneous  solution.  Hence  this  singularity  appears 
when  the  homogeneous  solution  results  in  a  vertical  wavenumber  that  fits 
perfectly  between  boundaries.  Yih  [36]  also  found  this  difficulty,  and  reports 
that  this  singularity  will  not  appear  if  the  analysis  is  extended  to  the  next 
order,  something  that  is  very  difficult  to  achieve,  and  is  not  pursued  here. 
At  this  time,  the  value  of  C2  near  7  =  invr  is  unclear. 

The  internal  waves  discussed  here  generate  a  mean  flow,  as  is  well-known. 
This  mean  flow  can  be  obtained  by  assembling  an  expression  for  S  using 
(133)  and  the  hnal  result  for  02,  then  averaging  over  a  horizontal  wavelength 
to  obtain  S.  The  hnal  mean  how  is  obtained  when  uq  is  combined  with  the 
vertical  derivative  of  S.  The  inhuence  of  this  wave-induced  mean  how  on  C2 
is  measured  by  the  quantity  T2  in  (152),  and  accounts  for  most  of  C2.  This 
was  also  the  case  for  Yih  [36].  Physically,  this  means  that  displacement  of 
streamlines  is  the  predominant  nonlinear  ehect,  rather  than  an  interaction 
of  higher  harmonics. 

Note  that  there  is  some  ambiguity  between  the  wavespeed  and  the  wave- 
induced  mean  how.  The  initial  steps  of  the  analysis  chose  a  coordinate  system 
moving  with  the  nominal  wavespeed  so  that  the  how  is  steady.  Any  other 
choice  of  coordinate  system  results  in  an  unsteady  velocity  held.  The  absolute 
speed  of  the  waves  can  only  be  dehned  relative  to  the  walls.  However,  as  the 
boundaries  have  no  shape,  and  the  horizontal  velocity  of  the  boundaries  is 
irrelevant  in  inviscid  how,  there  is  no  method  of  determining  the  absolute 
wavespeed.  The  second-order  correction  to  Cq  is  given  by  Q^C2.  Since 
is  positive  dehnite,  then  the  sign  of  this  correction  is  determined  by  the 
sign  of  C2.  Furthermore,  Q  contains  the  small  parameter,  e,  and  cannot  be 
larger  than  the  linear  wavespeed.  This  means  that  the  wave-speed  correction 
causes  the  wave  pattern  to  move  faster  or  slower,  depending  only  on  the 
wave  amplitude  and  the  sign  of  C2.  This  feature  still  represents  amplitude 
dispersion,  despite  this  ambiguity. 
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A  vertical  slice  of  6  at  second  order  for  an  example  set  of  parameters  shows 
that  this  wave  envelope  is  a  second-order  saturated  wave. 
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Figure  8:  Second-order  envelope  shape 


5  Nonlinear  Schrodinger  equations  for  a  wave 
packet  at  the  interface 

Consider  two  layers,  with  each  layer  distinguished  by  the  value  of  the  Brunt- 
Vaisala  frequency,  assumed  constant  within  each  layer.  The  density  prohle 
is  continuous  at  the  interface.  The  flow  is  modeled  with  the  Boussinesq 
approximation. 

A  packet  of  internal  waves  approach  the  interface  between  the  two  layers 
from  below.  The  waves  are  horizontally  periodic.  The  weakly  nonlinear 
theory  is  considered  below. 


5.1  Governing  equations 

Assume  incompressible  flow,  and  neglect  any  diffusion.  The  flow  is  assumed 
to  be  incompressible,  inviscid,  and  two-dimensional.  Stratihcation  is  present 
due  to  the  presence  of  a  non-diffusing  quantity,  and  the  flow  is  assumed  to  be 
Boussinesq.  The  flow  is  then  governed  by  the  Euler  equations  in  Boussineq 
form,  the  continuity  equation,  and  the  equation  of  incompressibility: 


Du  dp 
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where  (u,  w)  are  the  velocities  in  the  {x,  z)  directions,  respectively,  p  is  the 
dynamic  pressure,  po  is  an  average  (constant)  density,  is  the  mean  den¬ 
sity,  p  is  the  density  perturbation,  and  g  is  the  gravitational  constant. 

It  is  convenient  to  put  the  governing  equations  in  a  more  convenient  form 
by  eliminating  all  variables  in  the  linear  terms  in  favor  of  w.  The  result  is 
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where 


Ai  =  UUx  +  WUz, 


(159) 
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^3  =  UW^  +  WW^, 


(160) 

(161) 


A4  =  UPa;  +  WPz, 

5.2  Interfacial  conditions 

The  mean  density  profile  is  chosen  to  be  continuous,  but  have  a  discontinous 
hrst  derivative,  such  that  the  fluid  exists  in  two  semi-inhnite  layers,  each 
layer  having  a  unique  value  of  the  Brunt- Vaisala  frequency.  The  interface 
between  the  two  layers  must  satisfy  interfacial  conditions.  There  are  two 
types  of  interfacial  conditions  on  a  material  line  separating  two  inviscid  layers 
of  fluid;  kinematic  and  dynamic.  The  kinematic  condition  in  an  inviscid  flow 
states  that  the  normal  velocity  of  the  material  line  is  equal  to  the  normal 
component  of  velocity  of  the  fluid.  The  dynamic  condition  states  that  the 
pressure  must  be  continuous  across  the  material  line. 

The  kinematic  conditions  are 

r]t  + u~rix  =  w~ ,  (162) 

rit  +  u^Vx  =  w^,  (163) 

which  hold  on  the  interface,  z  =  p,  where  u~,w~  are  velocities  in  the  lower 
layer,  u^,w~^  are  velocities  in  the  upper  layer,  and  p  is  the  vertical  displace¬ 
ment  of  the  interface. 

A  primary  difficulty  is  meeting  the  interfacial  conditions  on  the  actual 
interface,  z  =  p,  without  knowing  the  position  of  the  interface  beforehand. 
This  difficulty  is  treated  by  expanding  all  terms  in  a  Taylor  series  about 
the  mean  position  of  the  interface,  in  the  same  manner  usually  used  for  free 
surface  flow.  The  kinematic  conditions  become 
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+  Ujp  +  -Uj^p  +  ... 

Vx  = 

X 

W+  +  w+p  +  -W+p^  +  ... 

,  (165) 

where  the  coefficients  are  now  evaluated  at  the  mean  position  of  the  interface, 
z  =  0. 

The  dynamic  condition  is  continuity  of  total  pressure,  p.  Hence 

p~  =  p'^  (166) 
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on  z  =  7],  where  p~  and  are  the  pressures  in  the  lower  and  upper  layers, 
respectively.  Consider  the  incompressible  flow  model  for  now,  rather  than  the 
more  restrictive  Boussinesq  approximation.  The  total  pressure  is  segmented 
into  the  mean  and  fluctuating  parts: 


p  =  p  +  p. 


(167) 


Expanding  (166)  in  a  Taylor  series,  as  before,  gives 


1 

+ 

1 

5  r-  -  1 

+  ^Ap  ] 

1  d"^ 

0+203  +  +  + 

II 

o 

o 

II 

[■p+  p+j 

II 

o 

+ 

+ 

+ 

+ 

1  3^ 

z=0 

+  ■■■  = 

+  ■■(168) 


Several  terms  may  be  eliminated  immediately.  Pressure  in  the  absence  of 
motion  is  continuous,  which  gives  p~  =  p'^  at  z  =  0,  allowing  these  terms  to 
be  dropped. 

Further  simplification  is  obtained  using  hydrostatic  equilibrium. 


dp 

dz 


-P  9, 


(169) 


(170) 


where  p~  and  p"*"  are  the  mean  densities  in  the  two  layers.  The  density  profile 
is  assumed  continuous  at  the  tropopause  (although  not  smooth),  implying 
p~  =  p'^  at  z  =  0.  Equations  (169)  and  (170)  then  result  in 


dp  dp'^ 

dz  dz 


(171) 


on  =  0.  These  two  terms  may  be  dropped  from  (168). 

Furthermore,  the  second  derivative  of  p  is  related  to  the  Brunt- Vaisala 
frequency: 


Pzz  =  -pzQ  =  -p 


Pz 

L  P 


=  +pN\ 


(172) 
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which  holds  for  each  layer.  Higher  order  terms  result  can  be  treated  in  a 
similar  manner.  The  dynamic  interfacial  condition  may  now  be  written  as 


1 

1 

+ 

_ 1 

+ 

Pz  -Pt 

1 

^+2 

Pzz  -  Ptz 

L  J 

2=0 

2=0 

+  ^Po 


N-^  -N+'^ 


2  I  1  PO 

V  +^  — 
3!  g 


N-^  - 


rj  -\ - 

r/3  +  ...  =  0,  (173) 


where  po  is  mean  density  at  the  interface. 

The  restriction  to  Boussinesq  flow  merely  implies  that  p  is  now  taken  as 
the  constant,  po. 


5.3  Weakly  nonlinear  theory 

Define  the  following  variables: 


—  X  Cpt^ 

(174) 

c 

=  ez, 

(175) 

r 

= 

(176) 

where  Cp  is  a  constant  to  be  determined  later,  and  e  is  a  small  parameter 
that  measures  the  vertical  length  of  the  wave  envelope.  Assume  all  variables 


depend  now  on  r,  and  2:. 

The  governing  equations  for  each  layer  become 

Po  [(m  —  Cp)  +  eUr  +  wUz  +  ewu;^]  =  — p^,  (177) 

Po  [{u  -  Cp)  +  eWr  +  wWz  +  cwwq]  =  -[Pz  +  epc]  -  PP,  (178) 

+  Wz  +  ewc_  =  0,  (179) 

{u  —  Cp)  +  epr  +  wpz  +  ewp(^  +  pzW  =  0,  (180) 

These  may  be  rearranged  into 

-poCpU^  +  epoUr  +p^  =  -PoFh,  (181) 

-poCpW^  +  epoWr  +pz  +  epc  +  pg  =  -poFy  (182) 

-Cpp^  +  epr  +  PzW  =  G.  (183) 
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Introducing  (174)  into  the  kinematic  interfacial  conditions,  (164)  and  (165), 
gives 


CpT]^  -  erjr 


u  +  [u^  +eu^  \7]  +  -  \u^^  +  eu^f.  +  ]rf  + 


w  +  [w^  +ew^  ]ri+  -  \w^^  +  + 


CpT]^  -  erfr 


u'  +  I  <  +  jr/  +  -  I  M+  +  eM+  +  + 


+  iwj  +  ew^  \r]+  +  ewj  +  r/^  +  ■ 


The  dynamic  condition,  (173),  becomes 


p  —  p 

1 

+  2 


+ 


2  =  0 


Pz  -  P  J  +  e  Pc  ~  Pc 


J  2=0 


Pzz-PJz  +(^[Pzc-PJc  [Pcc-Pic 


p  + 


2=0 


+  2^0 


iV-2  -iv+' 


2  I  1  Po 

p  +^— 
3!  g 


N-^  -  N+* 


?7  H - = 


5.4  Expand 

Consider  the  Fourier  expansions  given  by 


P  = 


imk^ 


m=—oo 

CO 


= 


^m(F^,C)e 


imk^ 


where  m  is  an  integer,  and  k  is  the  horizontal  wavenumber. 
The  operator  on  the  left  hand  side  of  (158)  becomes 


r±  (  d  d 


2  r 


(  (92  ,92 


+  iv±^41. 


■■  Pc 

,  (184) 

■■  Pc 

,  (185) 

0.  (186) 

(187) 

(188) 

(189) 
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Using  (188)  in  (158)  and  separating  Fourier  components  results  in 

(190) 

The  interfacial  boundary  conditions  are 

CpVm^  -  erimr  +  (191) 

Pm  Pm  ~  Qmi  (192) 

on  =  0. 


5.5  Linear  solution 


An  incident  wave  is  assumed  to  propagate  with  upwards  group  velocity  in 
the  lower  layer,  impinging  upon  the  interface.  The  vertical  component  of 
velocity  for  such  a  wave  is 


=  aA  +  cc, 


(193) 


where  A~  =  A~{(,t)  is  the  incident  wave  amplitude,  k  and  n~  are  the 
horizontal  and  vertical  wavenumbers,  respectively,  and  cc  means  complex 
conjugate.  The  resulting  dispersion  relation  is 


c^  = 


N 


_2 


(194) 


^  +  n~‘^ 

However,  (193)  by  itself  cannot  meet  the  interfacial  conditions.  A  reflected 
wave  as  well  as  a  transmitted  wave  is  required.  The  combination  of  waves  is 
given  by 

wi  =  +  cc,  (195) 

wf  =  +  cc,  (196) 

where  B~  =  B~{(,t)  is  the  reflected  wave  amplitude,  B~^  =  B^{(,t)  is 
the  transmitted  wave  amplitude,  is  the  vertical  wavenumber  in  the  upper 
layer,  and 

n~  —  vA 


K  = 

J  = 


n~  + 
2n~ 

n~  +  n+ 
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The  corresponding  displacement  of  the  interface  is 


r]i  =  ia- — J 
kc„ 


(197) 


where  A  is  the  amplitude  of  the  interface. 

5.6  Amplitude  equation 

The  left-hand-side  of  (190)  is 


^de  \de  dzy  de 

^  f  d^\  ^2 

I  7T7TT  +  TTW  I  +  2c„ 


Ord^yd^'^  dz"^  J  '^d^‘^dzdC,\ 


+  0(e^).  (198) 


If  Cp  is  chosen  to  be 


c^  = 


N 


±2 


^  ’ 

and  the  linear  solution  is  inserted  into  (190),  then  the  zeroeth  order  terms 
in  (198)  are  exactly  zero,  leaving 


^  92/92  92  ,  ^ 

—  2c„t^^:7T  I  TTW  +  TTW  +  2c: 


9r9/  V  '  9z'^ 
The  result  for  the  upper  layer  is 


94 


^deOzdC 


+  0(e^). 


aei 


kcp  ( /c2  n+  )  5+  -h  n'^k^clB^ 


2J^m-n+z)  ^ 


to  order  e2Q;.  For  the  lower  layer. 


aei 


kcj,  {k'^  +  n  +  n  k'^CpA^ 


2gi(fc«-n  2) 


aei 


kc„  (  A;2  +  n-^ )  B-  -  n-k^clB-  +  cc.  (199) 


The  linear  terms  in  the  interfacial  conditions  with  m  =  1  are  0('^)) 
do  not  add  to  zero,  but  instead  provide  a  relationship  between  the  incident 
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wave  amplitude  and  reflected  and  transmitted  wave  amplitudes.  The  nonlin¬ 
ear  terms  in  the  interfacial  conditions  are  no  choice  of  e  makes 

the  linear  and  nonlinear  terms  balance.  Hence  the  nonlinear  terms  in  the 
interfacial  conditions  with  m  =  1  are  higher  order,  and  the  linear  conditions 
apply: 


A- 

=  B- 

= 

C=o 

C=o 

(200) 


The  nonlinear  terms  in  the  interfacial  conditions  with  m  =  0  and  m  =  2 
both  contribute  through  wf  and  Mq  . 

The  nonlinear  terms  result  in 


=  0{ea^). 


The  nonlinear  terms  in  the  interfacial  conditions  result  in 

1  + 


n  n 

H  =  -2a^  —  J 


kCr, 


_2 


n 


Q  =  -4a  Pq-^K 


The  components  for  m  =  0,  2  can  be  extracted  by  inspection: 

Ho  = 


H2  =  -2a^^J 

kc„ 


T) 

Qo  =  Sa'^po—j^KAA*  O(ea^). 


_2 


Q2  =  -da  po 


2. 


k^ 


(201) 

(202) 

(203) 

(204) 

(205) 

(206) 
(207) 


5.6.1  Second  harmonic 

Equation  (190)  with  m  =  2  is  now 

L^(w2^=Mi,  (208) 

with  boundary  conditions 

CpP2^  -  ep2r  +  W2  =  H^,  (209) 
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pj  -ft  =  ft. 


(210) 


on  =  0. 

As  the  governing  equation  is  homogeneous,  the  only  particular  solution  is 
associated  with  the  boundary  conditions,  and  is 

^_2  r 

j^2 ^i(2k^^n2  z)  _j_  j^=^2 ^—i(^2k^zfn2  z) 


8n 


Wr,  =  a 


kcp{n2  +  n^) 


K 


where 


nt  =  n*-^  -  3k^. 


(211) 

(212) 


5.6.2  Mean  flow 


The  wave-induced  mean  flow  is  determined  using  the  wave  action  principles, 
as  outlined  by  Acheson  [1],  The  wave  energy  is 


E 


u'^  +  w'^  +  N'^rf 


(213) 


where  rj  is  the  displacement  of  a  material  line.  Again,  using  the  linear  solu¬ 
tions. 


E  =  2a^po 


n  ‘^  +  k‘^ 

P 


A-A-*  +  K^B-B- 
K{  A-*B-E‘^^~^  +  A-B-*e-^‘^^~^ 


E~^  =  2q;^po 


-I- 


k^ 


j2b+B+* 


The  mean  flow  is  determined  using 


E 


The  result  is 


Uq  =  2a  — 


-p  L 


n  +  k"^ 

P 


Uq  =  2q;  — 


Po'Wo  —  — • 

Cp 

A-A-*  +  K^B-B-* 

K[  A-*B-E‘^^~^  +  A-B-*e-^‘^^~^ 
j2b+B  +  * 


+  _o.2l  /n+VP 


/fc2 


(214) 

(215) 


(216) 


(217) 

(218) 
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5.6.3  Ml 

To  leading  order,  Mi  is 
Mf  =  -Aa^k^n-^  K‘^  5“  ^ 

+  K  A-  ,  (219) 


M+  =  0. 

The  result  is  three  amplitude  equations: 


P  +  iV 


=  0’ 


Cp(A;2  +  n-  ) 


/  X 

V  "  A;2  +  iV-2  W 


k“n  ^ 

- A-  B-  =  0, 

Cp(F  +  n-^) 


Recognizing  that 


=  0. 

+  n+^  ^ 


Then  these  amplitude  equations  are  more  generally  written  as 


9  d- 


i[A- -CgA-^] B-  A-  =  0, 


(220) 

(221) 

(222) 

(223) 

(224) 


/  \  71  C  ^ 

i(B;  -CgB^]  +A^K  A-  B-=  0,  (225) 

R+  +  CgB^  =  0.  (226) 

Note  that  Cg  is  the  vertical  group  velocity,  while  Cp  is  the  horizontal  phase 
velocity. 
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5.7  Results 


As  the  nonlinear  part  of  (226)  is  exactly  zero,  then  (226)  is  independent, 
linked  to  the  solution  of  (224)  and  225)  only  throught  the  interfacial  con¬ 
ditions.  Thus  (224)- (225)  may  be  treated  simultaneously  and  waves  in  the 
upper  layer  determined  separately,  if  needed.  It  is  also  possible  to  show 
several  constants  throughout  the  motion,  for  example. 


d 


(227) 


and  other  results. 

It  is  convenient  to  treat  (224)- (225)  numerically  using  the  leap-frog  method. 
The  results  show  that  the  interaction  of  the  incident  wave  with  its  relection 
beneath  the  interface  generates  a  higher  wave  amplitude,  and  corresponding 
mean  flow,  uq.  This  mean  flow  does  not  appear  above  the  interface,  and  the 
mean  flow  is  in  fact  discontinuous  at  the  interface.  Further  results  will  be 
submitted  soon  to  The  Journal  of  Fluid  Mechanics. 
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6  Stability  with  a  jump  in  stability  and  ve¬ 
locity 

This  problem  is  concerned  the  dynamics  of  waves  interacting  with  an  inter¬ 
face  that  has  both  a  jump  in  velocity  as  well  as  a  jump  in  the  Brunt- Vaisala 
frequency.  This  problem  is  motivated  by  recent  observations  ^  of  large  am¬ 
plitude  oscillations  over  Hawaii.  The  observations  remain  unexplained,  but 
suggest  resonant  over-reflection  of  internal  waves.  Previous  theoretical  results 
of  resonant  over-reflection  indicate  a  range  of  wavenumbers  where  waves  will 
spontaneously  be  created  by  a  shear  flow;  no  incident  wave.  These  results 
have  appear  in  an  AIAA  conference  paper  and  are  the  subject  of  further 
work  by  a  Ph.D.  student. 


6.1  Governing  equations 

Assume  incompressible  flow,  and  neglect  any  diffusion.  The  flow  is  assumed 
to  be  incompressible,  inviscid,  and  two-dimensional.  Stratihcation  is  present 
due  to  the  presence  of  a  non-diffusing  quantity,  and  the  flow  is  assumed  to 
be  Boussinesq.  Without  loss  of  generality,  a  coordinate  system  is  chosen  to 
be  moving  with  the  average  velocity  of  the  two  layers.  In  this  moving  but 
inertial  coordinate  system,  the  speed  of  the  upper  layer  is  U,  and  the  lower 
layer  is  —U,  following  Grimshaw  [JFM,  1979].  The  flow  is  then  governed 
by  the  Euler  equations  in  Boussinesq  form,  the  continuity  equation,  and  the 
equation  of  incompressibility: 


Po 


du 


^dit 


w 


du 

dz 


(228) 


Po 


dw 


dp 


^^dw  ^dw  dw 
±  G—  +u—  +  w— 
dx  dx  dz 


dp 

=  -¥z-P^^ 


,  ^^dp  ^dp  dp  dp 

±  U  —  h  u—  h  w—  h  —  0. 
dx  dx  dz  dz 


du  dw 
dx  dz 


0, 


(229) 

(230) 

(231) 


■^McHugh,  J.  P.,  Dors,  L,  Jumper,  G.  Y.,  Roadcap,  J.  R.,  Murphy,  E.  A.,  and  Hahn, 
D.  C.,  to  appear  in  JGR,  2008. 

^McHugh,  J.  P.,  Paper  AIAA-2009-0109,  47th  AIAA  Aerospace  Sciences  Meeting,  Or¬ 
lando,  FL,  January,  2009. 
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where  (m,  w)  are  the  velocities  in  the  {x,  z)  directions,  respectively,  p  is  the 
dynamic  pressure,  po  is  an  average  (constant)  density,  'p{z)  is  the  mean  den¬ 
sity,  p  is  the  density  perturbation,  and  g  is  the  gravitational  constant. 


6.2  Interfacial  conditions 


The  mean  density  prohle  is  chosen  to  be  continuous,  but  have  a  discontinous 
hrst  derivative,  such  that  the  fluid  exists  in  two  semi-inflnite  layers,  each 
layer  having  a  unique  value  of  the  Brunt- Vaisala  frequency.  The  interface 
between  the  two  layers  must  satisfy  interfacial  conditions.  There  are  two 
types  of  interfacial  conditions  on  a  material  line  separating  two  inviscid  layers 
of  fluid;  kinematic  and  dynamic.  The  kinematic  condition  in  an  inviscid  flow 
states  that  the  normal  velocity  of  the  material  line  is  equal  to  the  normal 
component  of  velocity  of  the  fluid.  The  dynamic  condition  states  that  the 
pressure  must  be  continuous  across  the  material  line. 

The  kinematic  conditions  are 


'nt-Ur]x  +  UiPa:  =  Wi,  (232) 

rit  +  Uria:  +  U2r]x  =  ^2,  (233) 

which  hold  on  the  interface,  z  =  p,  where  ui,wi  are  velocities  in  the  lower 
layer,  U2,W2  are  velocities  in  the  upper  layer,  and  p  is  the  vertical  displace¬ 
ment  of  the  interface. 

A  primary  difflculty  is  meeting  the  interfacial  conditions  on  the  actual 
interface,  z  =  p,  without  knowing  the  position  of  the  interface  beforehand. 
This  difficulty  is  treated  by  expanding  all  terms  in  a  Taylor  series  about 
the  mean  position  of  the  interface,  in  the  same  manner  usually  used  for  free 
surface  flow.  The  kinematic  conditions  become 


Pt  -  Upx 


ui+ui,p+  2'^u.V  + 


Wi  +  wi^p  +  -wu,p  + 


,  (234) 


Vt  +  Upx  + 


U2  +  U2,P  +  + 


W2  +  W2,p  +  -W2,,p  + 


,  (235) 
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where  the  coefficients  are  now  evaluated  at  the  mean  position  of  the  interface, 

z  =  0. 

The  dynamic  condition  is  continuity  of  total  pressure,  p.  Hence 

Pi  =  P2  (236) 

on  z  =  f],  where  pi  and  p2  are  the  pressures  in  the  lower  and  upper  layers, 
respectively.  Consider  the  incompressible  flow  model  for  now,  rather  than  the 
more  restrictive  Boussinesq  approximation.  The  total  pressure  is  segmented 
into  the  mean  and  fluctuating  parts: 


p  =  p  +  p. 


(237) 


Expanding  (236)  in  a  Taylor  series,  as  before,  gives 


[Pi  +  Pi] 


[p2  +  p2. 


z=0 


<9  1 

+  5zlp.+p.l 

a  ,,  . , 


z=0 


1 

2  dz‘^ 


V  +  ^^[Pi+Pi]  h^  +  ---  = 


z=0 


z=0 


1  <92 

P+2^1P2+P2] 


2  =  0 


r/2  +  ...  .(238) 


2=0 


Several  terms  may  be  eliminated  immediately.  Pressure  in  the  absence  of 
motion  is  continuous,  which  gives  pi  =  p2  at  z  =  0,  allowing  these  terms  to 
be  dropped. 

Further  simplihcation  is  obtained  using  hydrostatic  equilibrium. 


dpi 

dz 


-Pi9, 


(239) 


dp2 

dz 


=  -P29, 


(240) 


where  pi  and  p2  are  the  mean  densities  in  the  two  layers.  The  density  prohle 
is  assumed  continuous  at  the  tropopause  (although  not  smooth),  implying 
Pi  =  P2  at  z  =  0.  Equations  (239)  and  (240)  then  result  in 


dpi  _  dp2 
dz  dz 


(241) 


on  =  0.  These  two  terms  may  be  dropped  from  (238). 


58 


Furthermore,  the  second  derivative  of  p  is  related  to  the  Brunt- Vaisala 
frequency: 

pzz  =  -pzQ  =  -p  =  +pN‘^,  (242) 

L  p\ 

which  holds  for  each  layer.  Higher  order  terms  result  can  be  treated  in  a 
similar  manner.  The  dynamic  interfacial  condition  may  now  be  written  as 


1 

to 

_ 1 

+ 

1 

to 

_ 1 

1 

^+2 

plzz  -  p2zz 

L  J 

z=0 

L  J 

z=0 

+  ipo  -N2^  -  ^2"  7  +  . . .  =  0,  (243) 

where  po  is  mean  density  at  the  interface. 

The  restriction  to  Boussinesq  flow  merely  implies  that  p  is  now  taken  as 
the  constant,  po. 

6.3  Amplitude  equation 

Dehne  the  following  variables: 


e  = 

PC  Cpt  ^ 

(244) 

X  = 

e{x  -  Cgt) , 

(245) 

C  = 

(246) 

r  = 

e^t, 

(247) 

where  Cp  and  Cg  are  constants  to  be  determined  later,  and  e  is  a  small  pa¬ 
rameter.  Assume  all  variables  depend  now  on  P,,  x,  C)  and  2:. 

The  governing  equations  become 

Po  [{±U  -  Cp)  it^  +  e  {±U  -  Cg)  7  +  -|-  euu^  +  e^Ur  +  wuz  +  e'^wu,^] 

=  -[Pi  +  ePx]  >  (248) 

Po  [{±U  -  Cp)  w^  +  e  {±U  -  Cg)  +  uw^  -|-  euw^  +  e^Wr  +  wwz  +  e^ww(^ 

=  -[pz  +  -  pp,  (249) 
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(±[/  -  Cp)  p^  +  e  {±U  -  Cg)  +  up^  +  eupy,  +  e^Pr  +  wp^  +  +  p^w  =  0, 

(250) 

+  euy.  +  Wz  +  =  0,  (251) 

where  the  positive  sign  is  used  for  the  upper  layer,  and  the  negative  sign  for 
the  lower  layer. 

Introducing  (244)  into  the  kinematic  interfacial  conditions,  (234)  and  (235), 
gives 

-  -e(^U  +  Cg^Px  + 

Vi  +  ehx 

•  ,  (252) 


+ 


Ml  +  (  J ^  +  2  ]v^  +  ■■■ 

Mil  +  I  Mii^  +  ^  ^  i^izz  +  +  e'^Wic^c^  )  p^  + 


{u  -  Cp)p^  +  e(u  -  Cg^Px  +  e^Pr 


+ 


M2  +  (  U2z  +  ^  2  (  +  ■  ■  ■ 


Mi2  +  (  m;2z  +  e^t(;2c  )  ^  +  2  (  ^^2^  +  e^M;2zc  +  ]v‘^  ^ - 


Vi  +  ehx 
,  (253) 


The  dynamic  condition,  (243),  becomes 
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to 
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+ 

L  J 

2=0 

Piz  -V2z]+e  Pic  -  P2i 


J  2=0 


+  \po 


plzz  -  p2zz  )  +  (  P12C  -  p2zi  )  +  I  Picc  -  p2a 


p"^  + 


J  2=0 


Ni^  -  N2^ 


2  I  1  Po 

V  +  ^  — 
3!  g 


NC  - 


+  ■■■  =  0. 


(254) 


6.3.1  Reduced  equations 

The  governing  equations  at  each  order  have  leading  terms  that  allow  the 
reduction  of  the  equations  to  one  equation  operating  on  the  vertical  velocity. 
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The  generic  form  is 


Po(yTU  +  PoQi,  (255) 

poiyTU  +  Cj}jw^-P:,-  pg  =  poQs,  (256) 

PoiyTU  +  Cp  j  -  p^w  =  Qi,  (257) 

+  =  (258) 

where  the  Q’s  are  the  collection  of  terms  of  lower  order  that  appear  in 
each  equation,  different  for  each  order.  Eliminate  p,  u,  and  p  by  cross- 
differentiated  (255)  and  (256),  then  substituting  with  (257)  and  (258).  The 
result  is 

(  T  7/  +  Cp)  \^w  -  ~^w  =  (  ^  7/  +  Cp)  [gs,  -  Qi. 

H - Qi  -f  r  i  7/  -f  Cp'j  Qbz-  (259) 

Po  ^  ^ 

Defining 

N'^  =  _^,  (260) 

Po 

Then  (259)  becomes 

^  T  77  -f  Cp^  -f  N'^w  =  ^  =F  77  -|-  CpJ  —  Q\^ 

H - g4  +  f  T  77  -|-  Cp\  Q'oz-  (261) 

Po  ^  ^ 


6.3.2  Linear  stability 

The  most  general  form  for  a  disturbance  that  can  meet  the  interfacial  con¬ 
ditions  is 

+  (262) 

VJ2X  = 

+  (263) 
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where  the  terms  in  (262)  containing  Ai  and  Bi  are  incident  upon  the  interface 
from  z  =  —oo,  while  the  terms  in  (263)  containing  C2  and  D2  are  incident 
from  2:  =  +CX). 

If  the  incident  waves  in  each  layer  are  supressed,  then  this  linear  solution 


is 

wn  =  Cl  (e,  C,  r)  +  Di  (,  r)  ,  (264) 

W21  =  ^2  (e,  C,  r)  +  B2  (e,  C,  r)  .  (265) 

The  interfacial  conditions  result  in 


(^U  —  Cp^Ci  +  +  Cp^A2  —  0, 


(^U  —  Cp^Di  +  (^U  +  Cp^B2  —  0, 

rii  {u  +  Cp j Cl  -n2{u  -  Cp^  A2  =  0, 

ni  {u  +  Cp  j  Di  -  n2(u  -  B2  =  0. 
In  matrix  form,  these  are 


■  (G  +  Cp) 

0 

{U  -  Cp) 

0 

/7l2\ 

0  {U  +  Cp) 

0 

(U  -  Cp) 

B2 

n2{U  +  Cp) 

0 

-ni{U  -  Cp) 

0 

Cl 

0  n2{U  +  Cp) 

0 

-ni{U  -  Cp)_ 

\DiJ 

Setting  the  determinant  of  this  matrix  to  zero  results  in 

r  1 2 


ni(U  +  Cp)^  +  n2(U  —  Cpf^ 


=  0, 


(266) 

(267) 

(268) 
(269) 


0. 


(270) 


Implying  that  the  quantity  inside  square  brackets  is  zero.  The  resulting  equa¬ 
tion  is  identical  to  equation  (1.10)  of  Grimshaw  [JPM,  1979],  who  considered 
the  same  problem  with  constant  N  throughout. 

The  governing  equation  in  each  layer  determine  the  vertical  wave  number 
in  each  layer: 


n 


2 

1 


Nl 


(U  +  Cp)2 


m 


(271) 


n 


2 

2 


7V| 


(U  -  cp)2 


m 


(272) 
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Combining  (271),  (272),  and  (270)  gives 


{u  +  CpY 


(U 


-  m2  +  (7/  -  Cp) 


—  m2  =  0 


(273) 


There  are  apparently  two  modes  that  satisfy  (273),  as  outlined  by  Grimshaw 
[JFM,  1979].  One  mode  has  Cp  =  0. 

Rearrange  (273),  take  a  square  to  eliminate  the  square  root,  expand,  sim- 
plyify,  and  rearrange  to  obtain 


'IN^-NjUcpy 

_8  7/2^2 


niv2  +  iv2  1 

'INf-Nf 

4  7/2m2 

- 

_8  7/2m2 

0.  (274) 


There  are  three  solutions  to  this  third-order  algebraic  equation.  Note  that 
if  Ni  is  assumed  equal  to  N2,  then  the  results  of  Grimshaw  [JFM,  1979]  are 
recaptured. 

One  of  Grimshaw’s  modes  was  Cp  =  0;  this  is  no  longer  a  solution.  As 
a  result,  the  equation  does  not  reduce  to  a  second-order  equation,  and  the 
results  are  not  so  easily  obtained.  However  this  same  solution  exists  for  a 
nonzero  value  of  Cp,  and  corresponds  to  the  unstable  mode  for  the  original 
Helmholtz  instability  without  a  density  jump.  Note  that  with  a  density 
jump,  there  is  a  long-wave  instability,  whatever  the  velocity  difference.  This 
long-wave  instability  dissappears  with  the  density  jump. 


6.4  Results  from  linear  stability 

The  linear  results  are  shown  in  figures  9  through  11.  Figure  9  shows  all 
the  roots  of  (274)  for  ^  =  2,  indicating  all  three  solutions.  Note  that 
the  lower  branch  of  corresponds  to  a  zero  value  of  c*,  and  is  analogous 
to  Grimshaws  mode  with  c  =  0.  All  of  these  solutions  do  not  satisfy  the 
radiation  conditions.  Figure  10  shows  the  same  case,  but  only  shows  the 
modes  that  do  satisfy  the  radiation  conditions. 

Figure  11  gives  the  results  with  A^i  =  N2,  identical  to  the  constant  N 
case  of  Grimshaw  [12],  shown  here  for  comparison.  Grimshaw  showed  that 
there  are  two  modes,  one  with  c  =  0.  The  second  mode  can  be  seen  in 
figure  11.  For  ^  this  mode  does  not  satisfy  the  resonance  conditions. 


63 


kU/Ni 


Figure  9:  All  complex  wave  speed  with  ^  =  2 

and  for  ^  this  mode  is  unstable  (note  the  values  of  Cj).  Resonant 

over- reflection  for  this  mode  occurs  ioi  \  <  ^  <  ^. 

This  same  behavior  of  the  second  mode  exists  when  N  is  not  constant,  as 
shown  in  figure  10  for  ^  =  2.  There  is  a  critical  value  of  ^  of  approximately 
0.55  beyond  which  the  flow  is  unstable.  This  critical  value  depends  on  the 
ratio  Below  this  critical  value,  the  results  show  resonant  over- reflection. 
Note  in  hgure  10  That  the  lower  branch  of  this  second  mode  extends  to  ^ 
of  zero,  indicating  over-reflection  for  very  long  waves.  The  second  mode  for 
constant  N  does  not  satisfy  the  radiation  conditions  for  such  long  waves,  and 
does  not  exist.  However,  the  hrst  mode  does  not  exist  for  ^  =  2,  whereas 
this  hrst  mode  does  exist  for  all  long  waves  with  constant  N .  Overall,  the 
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Figure  10:  Complex  wave  speed  with  ^  =  2 


interval  of  wavenumbers  for  resonant  over-reflection  with  the  sudden  change 
in  N  is  somewhat  narrower  compared  to  the  case  with  constant  N.  Also  note 
that  only  modes  with  positive  values  of  Cr  exist  for  the  second  mode. 

A  case  where  N  decreases  with  the  vertical  is  shown  in  hgure  12  for  ^ 
=  For  this  case,  the  interval  for  resonant  overreflection  is  much  reduced, 
as  the  critical  value  of  ^  is  now  approximately  0.275.  Furthermore,  the 
wavespeeds  for  the  second  mode  are  all  negative,  opposite  the  case  with  N 
increasing  with  the  vertical. 
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7  Numerical  simulation  of  mountain  waves 
with  rotation 

The  third  continuing  task  is  treating  the  effects  of  the  Earth’s  rotation  on 
a  field  of  mountain  waves,  in  collaboration  with  T.  R.  Akylas.  The  work 
is  aimed  at  verihcation  of  a  recent  theoretical  result  of  Akylas  and  his  as¬ 
sociates.  The  computation  is  very  challenging,  as  the  physical  scale  of  the 
computation  must  be  very  large,  to  include  the  mountain  wave  effects  in  the 
far  held,  yet  must  also  accurately  determine  the  how  directly  over  the  moun¬ 
tain.  In  order  to  achieve  these  disparate  scales,  a  pseudo-three-dimensional 
approach  is  employed,  where  the  out-of-plane  velocity  is  non-zero,  but  does 
not  vary  in  this  third  direction.  Also,  the  domain  is  rectangular,  and  the 
mountain  is  included  with  an  artihcial  velocity  boundary  condition.  For 
mountains  of  small  height,  the  results  have  accurately  reproduced  Long’s 
solution.  However,  A  test  case  is  shown  in  hgure  (13). 

However,  as  the  mountain  height  is  increased  to  the  value  where  nonlinear 
ehects  are  important,  the  artihcial  boundary  condition  becomes  unstable.  A 
nonlinear  version  of  the  boundary  condition  has  been  implemented,  but  also 
is  unstable.  Promising  results  are  just  now  being  found  with  a  version  of  the 
immersed  boundary  method. 


8  A  vortex  pair  impinging  on  the  interface 

The  tropopause  is  associated  with  a  sudden  change  in  the  buoyancy  fre¬ 
quency,  which  is  a  quantity  that  is  intimately  related  to  internal  waves.  As  a 
result,  the  complicated  behavior  that  appears  near  the  tropopause  is  usually 
associated  with  internal  waves,  and  probably  usually  is  driven  by  internal 
waves.  However,  other  phenomena  can  be  effect  by  this  sudden  change  in  N, 
and  one  important  process  is  considered  here:  a  vortex  pair. 

The  tropopause  may  be  a  barrier  for  turbulence,  as  well  as  internal  waves, 
which  may  account  for  some  of  the  observations.  The  simplest  case  of  a 
turbulent-like  flow  is  a  vortex  pair.  The  vortex  pair,  or  vortex  ring  in  three 
dimensions,  will  propagate  with  a  fixed  speed,  raising  the  question  as  to 
whether  a  vortex  pair  will  reflect  off  of  the  tropopause  in  a  manner  similar 
to  internal  waves. 

To  investigate  this  possibility,  a  vortex  pair  has  been  studied  numerically. 
The  governing  equations  are  the  anelastic  equations.  The  numerical  method 
is  an  explicit  third-order  Adams-Bashforth  method,  spectral  in  space.  The 
work  was  performed  by  a  graduate  student,  Nicolas  Jenkins. 

8.1  Results 

A  vortex  pair  oriented  to  propagate  vertically  will  move  vertically  with  ap¬ 
proximately  constant  speed.  In  a  viscous  flow,  the  vortex  pair  gradually  loses 
energy,  but  maintains  it’s  coherent  form.  In  a  stratified  flow,  a  strong  vortex 
pair  will  act  mostly  like  the  constant-density  case,  however  a  weak  vortex  or 
vortex  pair  will  disintegrate  into  internal  waves  very  quickly.  In  a  viscous 
stratihed  flow,  an  initially  strong  vortex  pair  will  act  like  the  constant-density 
case  until  the  pair  has  lost  some  energy,  and  then  the  vortex  pair  will  disin¬ 
tegrate  into  internal  waves. 

The  research  considered  a  relatively  strong  vortex  pair  such  that  the  vortex 
pair  retains  a  coherent  form.  The  pair  was  released  below  the  interface  and 
allowed  to  impinge  upon  it.  The  results  are  shown  in  figure  (14)  for  three 
time  steps,  with  the  top  figure  the  earliest  time  and  the  bottom  the  latest. 
The  top  hgure  shows  the  vortex  pair  beneath  the  interface  (shown  with  the 
dashed  line),  the  middle  hgure  the  pair  has  passed  through  the  interface, 
and  in  the  third  hgure,  the  pair  has  returned  to  the  interface.  What  happens 
is  that  the  vortex  pair  is  trapped  by  the  interface  and  winds  up  oscillating 
through  the  interface.  Note  that  in  a  layer  of  constant  N,  the  vortex  pair 
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continuous  to  propaga  upward  until  internal  waves  are  created. 

These  results  show  that  indeed  the  tropopause  is  a  hlter  for  certain  scales  of 
vortical  motion,  not  just  internal  waves.  It  is  likely  therefore  that  turbulence 
is  hltered  by  the  tropopause,  which  may  account  for  some  of  the  observations. 
These  results  will  soon  be  submitted  to  The  Journal  of  Fluid  Mechanics. 
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